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ABSTRACT 


The spreading behavior of new or invasive species is a central topic in ecology. The 
modelings of free boundary problems are widely studied to better understand the 
nature of spreading behavior of new species. From mathematical modeling point of 
view, it is a challenge to perform numerical simulations of free boundary problems, 
due to the moving boundary, the stiffness of the system and topological changes. 

In this work, we design numerical methods to investigate the spreading behavior 
of new species for a diffusive logistic model with a free boundary and a diffusive 
competition system with free boundaries. We develop a front-tracking method, which 
explicitly tracks the location of the moving boundary, in one dimension and higher 
dimensions with spherical symmetry. In higher dimensional cases, we introduce level 
set method to handle topological bifurcations. Various numerical simulations in one 
and two dimensional spaces are presented to validate the accuracy, and stability of the 
proposed numerical methods. To efficiently solve stiff reaction-diffusion equations, we 
also develop implicit integration factor (IIF) method combining Krylov subspace to 
solve the diffusive logistic model with a free boundary in one dimension. Compared 
with different numerical schemes, it can be observed that Krylov IIF is advantageous 
to other approaches in terms of stability and efficiency by direct comparison through 


numerical examples. 
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CHAPTER 1 


INTRODUCTION 


The spreading behavior of new or invasive species is a central topic in ecology, and 
considerable research has been devoted to the better understanding of the nature 
of such spreading. The modeling of biological invasions has been widely studied in 
[1, 2, 19, 33, 37, 54, 61, 62] under the crucial restriction that the spatial domain is 
not constrained by the population behavior. The first diffusive logistic model related 
to biological invasions was initiated in 1937, of course without boundary restrictions, 
independently by Fisher [19] and Kolmogorov-Petrovsky-Piskunov (KPP){33]. To 
our knowledge the seminal paper [14] by Du and Lin is the first contribution in the 
field of spreading of populations where a Stefan condition is used and managing a 
moving boundary problem of parabolic type. Further developments of this problem 
have been treated in [13, 15, 16]. The diffusive logistic model with a free boundary 
of [14] for the density of population of the invasive species U(t,x) depending on time 


t and spatial variable x states as follows: 


“ — DAU = f(U) for x E€ M(t), t>0; U=0 forx €ON(t),t>0. (1.1) 


The nonlinear function f(U) is assumed to be a C! function satisfying f(0) = 0, 
and in the literature it is often taken to be the logistic function f(U) = U(a — bU) 
with a,b positive constants. In our work, we take the logistic function as f(U) to 


demonstrate the numerical methods. 


The evolution of the moving domain (t) C RY, or rather its moving boundary 


O§2(t) is determined by the one phase Stefan condition, which, in the case 0{2(t) is a 


C' manifold in R”, can be described as follows: 


Any point x € 02(t) moves with velocity pu|V,U(t, x)|n(x), where n(x) 


is the unit outward normal of (2(t) at x, and p is a given positive constant. 


The moving boundary O{2(t) is generally called the “free boundary”. The spreading 
behavior in the diffusive logistic model with a free boundary has been extensively stud- 
ies theoretically. Specially, the authors in [14] establish the spreading-vanishing di- 
chotomy in one dimension: Suppose (u(t, x), h(t)) is the solution of the free boundary 
problem in one dimension, either new species spreads successfully so that h, = +0o 
and 


lim (ult, i) uniformly for x in any bounded set of [0, 00) 


or new species vanishes as 


7 1D 
ho Soi 7 and lim, |lu(t, )Iccon@y = 9. 


Theoretical results in [14] also state that for hy > 2/2 spreading of the species 
is guaranteed. Even if ho < 2/2, spreading occurs under condition pp > ~* where 
w* is an unknown threshold depending on u(0,2), see Theorem 3.9 of [14]. In the 
spreading case the population density tends to the habitat carrying capacity limit ¢ 
as time tends to infinity, see Lemma 3.2 of [14]. For ho < a2 and js < y* vanishing 
happens, satisfying that z,/ = is an upper bound of h(t). 

In [16], the regularity and long-time behavior of O0(2(t) and u(t,x) in high di- 


mension are studied. It is shown that a spreading-vanishing dichotomy holds: either 


Q(t) stays bounded (i.e., is contained in some fixed ball in R%) for all t > 0, and 


u(t,x) > 0 as t > oo uniformly in x € Q(t), or Q(t) converges to RX, with OQ(t) 


approximating a moving sphere enlarging to infinity as t + oo. Moreover, in the 
latter case, for all large t, O(2(t) is a smooth closed manifold without boundary. 

In contrast, very few papers have treated numerically these nonlinear models 
focusing on the stability and the preservation of the qualitative properties of the 


theoretical solution [14, 16]. On one hand, extremely small time steps generally are 


required due to the stiffness of the system. On the other hand, it is always difficult 
to efficiently and accurately handle the moving boundaries. 

To efficiently handle the moving boundaries, level set method [18, 45, 51, 52, 65, 66] 
and front tracking method [35, 47, 63, 67] are two popular numerical approaches. One 
distinct feature of front tracking [11, 22, 27, 30, 36, 58] is using a pure Lagrangian 
approach to explicitly track locations of interfaces, but it is difficult to handle topo- 
logical bifurcations in high dimensions, while the level set method can efficiently 
overcome such difficulties. 

In this work, we are devoted to designing efficient numerical methods to solve the 
diffusive logistic model with a free boundary. Combining the implicit solver, front 
tracking method, which explicitly tracks the location of interfaces, is adopted to 
develop numerical schemes for the one dimensional case and higher dimensional case 
in a spherically symmetric setting. This approach enjoys the following advantages: 
(1) overcome the stiffness of the system caused by diffusion term and reaction term, 
so large scale time steps are allowed; (2) capture the locations of the moving front 
explicitly. A front-fixing approach is also introduced and applied to develop numerical 
schemes for the one dimensional case and higher dimensional case in a spherically 
symmetric setting to verify the accuracy of front-tracking method. The main idea 
of front-fixing approach is to transform the original moving boundary problem into 
a problem with a fixed computational domain. Then a new but equivalent system 
is obtained. Comparison for front-tracking method and front-fixing method are hold 
in numerical experiments to show accuracy and consistency of both methods. In 
higher dimensional case, we introduce level set method to handle moving boundaries 
instead. As front-tracking method tracks locations of interfaces explicitly, it is difficult 
to handle topological bifurcations in high dimensions, while level set method can 
efficiently overcome such difficulties. The idea behind level set method is to construct 


a level set function @ and move @¢ with the correct speed at the front and followed by 


updating u(t,x). The new position of the front is stored implicitly in ¢. With this 
approach, we avoid difficulties that arise from explicitly tracking the moving front and 
thus increase the efficiency to deal with complex interfacial geometries. Numerical 
experiments are designed to illustrate and verify the theoretical solution [14, 16]. 
We then extend our numerical methods to the diffusive competition system with 
two free boundaries to investigate the spreading behavior of two computing species. 
The rest of this dissertation is organized as follows. In Chapter 2, we briefly 
review the diffusive logistic model with a free boundary. A front-tracking approach 
is introduced for one-dimensional model and for a two-dimensional case with radial 
symmetry, and the method is also compared to front-fixing method. After that, a 
level set method is discussed for a more general two-dimensional case for one-species 
reaction-diffusion system with a free boundary. Finally we present numerical simu- 
lations to validate the accuracy, efficiency, stability and consistency of the proposed 
numerical methods for the diffusive logistic model with a free boundary. In Chapter 
3, we develop front tracking and front fixing method for the diffusive competition 
system with two free boundaries in one dimension, and we also introduce the level 
set method for the diffusive competition system with two free boundaries in the 
general two-dimensional case. The accuracy, efficiency, stability and consistency of 
the proposed numerical methods for the diffusive competition system with two free 
boundaries are tested. In Chapter 4, we first transform the one-dimensional diffusive 
logistic model with a free boundary into a system with a fixed computational domain, 
and then introduce four different temporal schemes: Runge-Kutta, Crank-Nicolson, 
implicit integration factor (IIF) and Krylov IIF for handling such stiff systems. Nu- 
merical examples are examined to illustrate the efficiency, accuracy and consistency 
for different approaches, and it can be shown that Krylov IIF is superior to other 
three approaches in terms of stability and efficiency by direct comparison. Finally 


some concluding remarks and future work are given in Chapter 5. 


CHAPTER 2 
NUMERICAL METHODS FOR A DIFFUSIVE LOGISTIC 


MODEL WITH A FREE BouNDARY! 


1g. Liu, Y. Du and X. Liu. Accepted by International Journal of Computer Mathematics. 
Reprinted here with permission of publisher, 11/13/2019. 


2.1 ABSTRACT 


The spatial-temporal spreading of invasive species in a habitat is a central topic in 
ecology. It is described by a diffusive logistic model with a free boundary, where the 
free boundary represents the unknown expending front of the species. We propose 
front-tracking method and level set method focusing on the preservation of the quali- 
tative properties of theoretical solutions. Illustration with numerical examples of the 
spreading-vanishing dichotomic behavior are given. The stability, accuracy and con- 


sistency of proposed numerical methods are established with numerical experiments. 


2.2. INTRODUCTION 


The spreading behavior of new species of the diffusive logistic model with a free 
boundary is extensively studied in recent years theoretically [13, 14, 15, 16]. In one 
dimension, the spreading-vanishing dichotomy in the diffusive logistic model with a 
free boundary establish that (see Theorem 3.3 in [14]): 

Let (u(t, xz), h(t)) be the solution of the free boundary problem in one dimension. 
Then the following alternative holds: 

Either 


e spreading: ho = +00 and limy,400(u(t,)) = ¢ uniformly for x in any bounded 


set of [0, 00); 
e vanishing: h, < x [2 and lims++00||u(t, Je qo,aw)) = 0. 


Theoretical results in [14] establish that for ho > 2/2 spreading of the species is 


guaranteed. Even if ho < 24/2 


x\ > Spreading occurs under condition > p* where 


u* is an unknown threshold depending on u(0,2), see Theorem 3.9 of [14]. In the 
spreading case the population density tends to the habitat carrying capacity limit ¢ 


as time tends to infinity, see Lemma 3.2 of [14]. 


In [16], the regularity and long-time behavior of 0Q(t) and u(t, x) in high dimen- 


sion are studied, and it is shown that a spreading-vanishing dichotomy holds: either 


Q(t) stays bounded (i.e., is contained in some fixed ball in R’%) for all t > 0, and 


u(t,x) > 0 as t > oo uniformly in x € Q(t), or Q(t) converges to RX, with OQ(t) 
approximating a moving sphere enlarging to infinity as t + oo. Moreover, in the 
latter case, for all large t, O(2(t) is a smooth closed manifold without boundary. 

We aim to treat the diffusive logistic model with a free boundary in one dimen- 
sion and higher dimensions numerically, which can be a continuation and numerical 
complement of [13, 14, 15, 16]. We propose front-tracking method in one dimension 
and higher dimension with radially symmetry, we also introduce level set method 
for general two dimension case. The dichotomic behaviors of new species from the 
diffusive logistic model with a free boundary are illustrated by numerical examples. 
We also test the stability, accuracy and consistency of proposed numerical methods 


with numerical experiments. 


2.3. NUMERICAL METHODS FOR 1D MODEL 


The diffusive logistic model with a free boundary in one dimension of [14] for the 
density of population of the invasive species U(t,x) depending on time t and spatial 
variable x states as follows: 


2 
OU p?2 U 


ae azz = Ula dv), (20,. Ose Ae), (2.1) 


together with the boundary conditions 
OU 


Dg Wt 9) =U, UG HG) =O. $30; (2:2) 
the Stefan condition 
; OU 
A(t) ==f—@, AO) tS 0; (2.3) 
Ox 
and the initial conditions 
HO) SHG <0, 6) = Oe). Oar Ag: (2.4) 


The initial function Up(x) satisfies the following properties: 
Uo(x) € C?([0, Hol), U3(0) = Uo( Ao) = (0, Uo(x) >0, O<a< Ao. (2.5) 


Here H(t) is the unknown moving boundary such that the population is distributed 
in the interval 0, H(t)]. D > 0 is the dispersal rate and the positive parameters a and 
6 are the intrinsic growth rate and the intraspecific competition, respectively. The 
parameter j > 0 involved in the Stefan condition (2.3) is the proportionality constant 


between the population gradient at the front and the speed of the moving boundary. 


2.3.1 A FRONT TRACKING METHOD FOR 1D MODEL 


In this approach, we combine the Lagrange method to track the interface with an 
implicit scheme to solve the parabolic equation to update U on all the grid points. 
Let us consider the moving front problem (2.1)-(2.4) in a fixed computational domain 
(0, 7] x [0, L], i-e., we want to learn about the distribution of population of the invasive 
species U(t,z) in the region [0, L] at time T. Set the spatial step size discretization 
h = Ax = L/M, and the time step size discretization k = At = 0.1Az. The mesh 
points (¢”,2,;), with t? = nk, n > 0, a; = ih, 0 <i < M, where M are the total 
number of subintervals of [0, Z]. Let us denote the approximate value of U(t", x;) at 
the mesh point (t",z;) as 


Ue UE a). (2.6) 
and let H” be the approximation of H(t"). 
Step 1. Track the position of the moving front. 


According to the Stefan condition (2.3) 


H'() = w(t. MO), t>0, 


Here we consider using the central approximation of the spatial derivatives to ap- 


proximate et, H(t)), which can be divided into the following four cases: 


1. When a; < A” < 2441, 1 = 2,3...M — 1, we denote d = ee Let us first 
consider the symmetric point of x;_; respect to the position H”, which is denoted by 
%;-1. Especially when H” = x;, %j;-1 = X11. We use the Lagrange extrapolation to 
construct polynomial P” from the value of d, h, u®_,, u?_,, u? and H”, thus at Zj_1, 


we use the value of P” at %;_; instead of u(t", %;_1), therefore 


OU PEGS) — y 1 

—(t", H”")& = p= 2,3,...,M@—-1. De 
aE ) ) 2(1+d)h y) a Ob ees d ( ) 
TQ Li—2 Lji-1 a lle Hiei. 2M 


(l+d)h  (1+d)h 


Figure 2.1: Illustration of how to evaluate 94(t, H(t)) at t” when x; < H” < wi4i, 
ie DAM Ri, 


Remark: The most challenging part of front tracking method is the evaluation of 
Mt, H(t)), where H(t) is moving along time. It causes the difficulty of finding a uni- 
form finite difference approximation to evaluate 24(¢, H(t)) with high order accuracy 
because the distances among the points x7;_;, x; and H” are disunion. The numerical 
methods of evaluating 5 8 (t, H(t)) on nonuniform points have some backwards, they 
usually cause singularity when the distances among H” and grid points x;, xj41 are 
quite small. The methods of evaluating 94(t, H(t)) by combining the evaluation of 
S(t, 2;) and 24(t, x41) can avoid the singularity, however, it disturbs the accuracy 
of the front tracking method when combined with the process of updating U(t, x). In 
our method, we choose the classical central approximation of the spatial derivatives 
combining the Lagrange extrapolation to evaluate 24(t, H(t)) on uniform points in- 
stead on nonuniform points to keep second order accuracy in space. And also, when 
x; < H” < 241, we evaluate a(t, H(t)) on z;_,; and Z;_; instead of x; and Z; to 
avoid singularity when H” is very close to x;. Actually, when H” = 2;, Equation 


(2.7) turns into 


which shows that our method has a good consistency. 

2. When 0 = x < H” < 21, the central scheme approximation of the spatial 
derivatives to approximate Mt, H(t)) involves the fictitious value u”, at the point 
(t”, —h). The value u”, can be estimated from the second-order discretization of the 


boundary condition (2.2), 


which implies that u", = uj = 0. It is obvious that all the values of u? on the 
grid points are equal to 0 except uj. Numerically, we take H'(t) = 0, and it can 
be explained that the species is only located inside one grid mesh. The simulation 
should stop here indicating that a more refined mesh is needed. 


3. When x, < H” < 22, we denote d = tL 


. Let us first consider the symmetric 
point of x9 respect to the position H”, which is denoted by %). Then we consider the 
value of u", = u'’, and use the Lagrange extrapolation to construct polynomial P” 
from the value of h, d, u",, u?, uj and H”. Then at Zo, we use the value of P” at % 


instead of u(t”, Zo). 


(1+ d)h (1+d)h 


Figure 2.2: Illustration of how to evaluate Mt, H(t)) at t” when x1 < H" < x9. 


4. When H” = xy, it implies that the spreading of the populations already 
goes out of the computational domain [0, Z], and the simulation should stop here 
indicating a larger computational domain is needed. 

Remark: In the system (2.1)-(2.4), case 2 and case 3 will never happen, since 
jt > 0, the front H(t) is increasing along time t. Actually, the front tracking method 


possesses preferable adaptability to track diversified conditions of front. 


Step 2. Update the value of U(t"'", x;). 
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1. When x; = H”*', we know that x; is the location of the boundary at ¢”*?. 
As the density of the population of new species at the boundary and beyond the 
boundary is 0, we set ut! = 0 and unt = 0, for j =i1+1,i+2,...M. Let us consider 
the central approximation of the spatial derivatives at x;, for 7 = 0,1,2,...,7 — 1, 


where U is updated by the backward Euler 


n+1 n n+1 n+l 
ug’ —U 20 20 
ane remes 2) 7B Oo + att) — b(untt)?. 
yeti — yr ynti Qyrtt as yrtl (2.8) 
j fod ar) j+1 1 1\2 . 
, i — p- 53 Pe taught — bury, 7 =1,..40-1. 


Then we use Picard Iteration (or Newton Iteration) to solve nonlinear system (2.8). 


2. When x; < H™*! < 2444, denoting R = yer at we use the Lagrange extrapolation 


to construct polynomial P’ from the value of h, R, uty, utt/, u?*t and H"*!. Then 


at Xi41, we use the value of P” at x4, instead of ui4', where U is updated by the 


backward Euler in time 


n+1 n n+1 n+1 
Ug Uo pa _ 2U9 n+1 b( aes 
: _ 72 | 0 U . 
n+l n n+1 n+1 n+1 
Tale Us — 2u;"° + Uy ’ 
= = pa fp OR OPE, FS 1nd = 1 (2.9) 
n+1 n n+1 n+1 L 
Ur — Uy wey — Qu") + PY (ai41) n+1 n+1)2 
i =D h2 tauy" — b(u;")°. 


Picard Iteration (or Newton Iteration) will be applied to solve nonlinear system (2.9). 


Remark: Front-fixing method by transformation [48] 


As discussed in [48], front-fixing method has been introduced to solve 1D diffusive 
logistic model (2.1)-(2.5) by transforming it into a problem with a fixed domain {0, 1]. 


Let us consider using the Landau transformation [12, 34], 


2 eS W(t,2) = UG,2). (2.10) 


a, 
A(t)’ 
Then the 1D diffusive logistic model (2.1)-(2.4) turns into the form: 


ow | 20W ew 


Gt) 
where G(t) = H?(t). 
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Boundary conditions (2.2) and Stefan condition (2.3) take the forms 


OW 


and 
Cw= =e ile ae ed (2:13) 


respectively, while the initial conditions (2.4) become: 
G(0) = H?, W(0,z) =Wo(z) =Uo(zHo), O<z< 1. (2.14) 
Conditions (2.5) for the initial function Up(x) are translated to Wo(z) as follows: 
Wo(z) € C?([0,1]), W5(0) =Wo(1) =0, Woz) >0, O<z<1. (2:15) 


After the transformation, the new problem lies in solving the nonlinear parabolic 
partial differential equations (2.1) in the fixed domain [0,7] x [0,1] for the variables 
(t, z). Let us consider the step size discretization k = At, h = Az = 1/M, and the 
mesh points (t", z;), with t’ =kn, n> 0, z; = jh, 0<j < M and M is the number 
of intervals in [0,1]. Let us denote the approximate value of W(t”, z;) at the mesh 
point (t”, z;), 


we = Wt", 2;), (2.16) 


and let g” be the approximation of G(t"). 


Considering the forward approximation of the time derivatives 


wrtt_w? aw ght! — g 
z ae) 2). = a ee), 217 
and the central approximation of the spatial derivatives, 
wr We OW on wha — 2WF + Wy | OW 
~ ieee & Ue) 2.18 
oh Oz ( 2%), h2 Oz2 ( 24) ( ) 


Equation (2.11) is approximated based on Equations (2.17) and (2.18), 


n+1 yn qr ccigit n+l n Te > 62 n n 
5 Wi 25 Witt wri 9 — 9 wir 205 Th 


k Z 2h 


g 


= gw; (a—bw;), n2=0,0<75M-1. 


Transformed Stefan condition (2.13) is discretized using first order forward approxi- 


mation for G’(t) and three points backward spatial approximation of oY tt, 1) 


; = —F(3why — 40h + wiya), n>0. (2.20) 


to preserve accuracy of order O(k) + O(h?). 


2.4 NUMERICAL METHODS FOR 2D MODEL WITH RADIAL SYMMETRY 


To solve the 2D diffusive logistic model in polar coordinates, the system can be written 
as 


au PU 10U ~, 1 8U 


Br = U(a-— <O< 
Ot (Fp r Or r2 79? U(a bU), t> 0, 0 S 7 _ 27, r> 0, 


where (rcos6@,rsin8) € Q(t). 

We assume that the environment and the solution are radially symmetric, i.e. we 
set the initial domain {2 is a disk, the initial function Up(x) is radially symmetric, 
the moving boundary O{2(t) is thus a circle whose radius we denote by H(t) and the 
solution U(t,r,@) = U(t,r), the 2D diffusive logistic model with radial symmetry can 


be written as a 1D diffusive logistic model 


OU 0’U 10U 


AE (apa? a an) — Ue OU), PSUs Uaeeore). (2.21) 


together with the boundary conditions 


OU 
3, (9) = (i, UG =0c oS, (2.22) 
the Stefan condition 
; OU 
A'(t) = ~hay(t, H(t), t>0, (2.23) 
and the initial conditions 
H(0)= Ho, U(0,r)=Ud(r), O<r< A(0). (2.24) 


13 


2.4.1 FRONT-TRACKING METHOD FOR THE 2D MODEL WITH RADIAL SYMMETRY 


We solve the 2D diffusive logistic model with radial symmetry (2.21)-(2.24) similarly 
to the 1D case. As the process shown above, we use an implicit scheme to solve the 
parabolic equation to update U on all the grid points. 

Let us consider the moving front problem (2.21)-(2.24) in a fixed computational 
domain [0,7] x [0, LZ], i.e., we want to find the distribution of the population in the 
region [0, ZL] up to time T. Set the step size discretization h = Ar = L/M,k = At = 
0.1Ar, and the mesh points (t”,r;), with t? = kn, n> 0, 7; = ih, O<i < M and 
M is the total number of subintervals of |0, Z]. Let us denote the approximate value 


of U(t",r;) at the mesh point (t”",r;) by 

ti PO rey (2.25) 
and let H” be the approximation of H(t"). 

Step 1. Track the position of the moving front. 


According to the Stefan condition (2.23), we use central approximation of the 
spatial derivatives to approximate at, H(t)): 

1. When r; < A” < rigi, 1 = 2,3,...,M —1, denoting d = aa firstly, we 
consider the symmetric point of r;_; with respect to the position H”, which is denoted 
by 7;-1. Specifically when H” = r;, 7;-1 = rj41. We use the Lagrange extrapolation 


method to construct a polynomial P” from the value of d, h, u% 5, u®4, u? and 


H™, and at #;_1, we use the value of P” at 7;_1 instead of u(t", 71), 


TO Misa Tisd Tm HH” i "yl 
TL a#.—-Y 


(lt+ ah (l4+dh 
Figure 2.3: Illustration of how to evaluate 2 (t, H(t) at ¢” when r; < A” < ryt, 


§=2,3..M—1. 


OU ae P¥(F_1) = Us 4 
~~ Dean «2 


10 3 MM, (2.26) 
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Remark: The most challenging part of the front tracking method is the evaluation 
of “t,H (t)), where H(t) is the moving boundary. It is difficult to find a uniform 
finite difference approximation of at, H(t)) with high order accuracy because the 
distance of the moving point H” to the set of grid points {r;} does not have a uni- 
form positive lower bound. The numerical methods of evaluating 24(t, H(t)) in such 
a case need to be carefully designed to avoid singularity caused by H” being very 
close to some grid point r;. Evaluating Mt, H(t)) by combining the evaluation of 
Mt, r;) and 2 “it, ri41) can avoid such singularity, however, it destroys the accuracy 
when combined with the process of updating U(t,r). In (2.26), we combine the clas- 
sical central approximation of the spatial derivatives and the Lagrange extrapolation 
method to evaluate Ht, H(t)) to ensure second order accuracy in space. For exam- 
ple, when r; < H” < rj41, we evaluate aU ~(t, H(t)) on r;_1 and F;_; instead of r; and 
7; to avoid singularity when H” is very close to r;, and when H” = r;, (2.26) becomes 


OU P¥(ri41) — ye 1 
aie. n H” ~ wm 
Or ye) 2h , 


(S03. 2 WR: 


2. When rp < H” < r1, the central approximation of the spatial derivatives to 
approximate Xt, H(t)) involves the fictitious value u", at the point (t”,—h). The 


value u", is eliminated from the discretization of the boundary condition (2.22), 


which means that u", = uf = 0, we can see that all the values of u’? on the grid 
points are equal to 0 except uj. The simulation should stop here indicating that a 


more refined mesh is needed. 


3. When r,; < H” < ro, denoting d = ~=—+, we consider the symmetric point of ro 
with respect to the position H”, which is denoted by 79. Then we consider the value 
of u", = u', and use the Lagrange extrapolation method to construct polynomial P” 
from the value of h, d, u",, u%, u} and H”, and at 7, we use the value of P” at 7% 


instead of u(t”, 70), 
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(1+ d)h (1+d)h 


Figure 2.4: Illustration of how to evaluate Ht, H(t)) at t” when r; < H" <1. 


4. When H” = ry, it means that the spreading of the populations goes out of 
the computational domain [0, LZ], and the simulation should stop here indicating for 
a larger computational domain. 

Remark: In the system (2.21)-(2.24), case 2 and case 3 will not happen, since 
the front H(t) is increasing in time t. However, front tracking method possesses 
preferable adaptability to track various moving front conditions, such as those used 
in [3], where the front need not be increasing in time, and therefore case 2 and case 


3 indeed occur. 
Step 2. Update the value of U(t"t', r;). 

1. When r; = H”*', we know that the boundary locates at r; at t’t!. Setting 
ut! — Oand Ge = 0, for 7 = i+1,i+2,...,.M, we consider the central approximation 


of the spatial derivatives at x,;, for 7 = 0,1,2,...,2 — 1, where U is updated by the 


backward Euler method in time 


n+l in n+1 n+1 
Uo Ug = pau 2U9 \ aunt — b(ugt?y? 
k h2 ! 0 0 : 
ntl 4,7 n+l n+1 n+l n+l n+l 
ee 2uj) + Ups | Uj+i MITT quttt — o(uttty? 
ro 2 P2 j j , 
k h 27h 


where j = 1,2,...,4 — 1. We use Picard Iteration (or Newton Iteration) to solve the 


nonlinear system. 


2. When r; < H"*! < r,4,, denoting R = aie we use the lagrange extrap- 


olation to construct polynomial P’ from the value of h, R, uty, uth, ult! and 


H™*' and thus at 7;,1;, we use the value of P” at r;,, instead of Was where U is 
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updated by the backward Euler method in time 


n+1 n n+1 n+1 
Up — Ug prt BMGT oe a Het — blurt) 
n+1 an n+l n+1 1 ,n+1 n+1 _ , n+l 
Up TF py Maan TU TT) Lmtd p(_yn thy? 
( y F ) r au, (u" ) : 
k h2 2jh2 j j 
n+l n n+1 nt+1 L L n+1 
UT D wey — 2ugh + PY (tin) | P' (ini) — Ua Laut! — b(yntty2 
k ca ( h2 Jih2 ) r QU; ~~ (u; ) : 


Where j = 1,2,...,4 — 1. Picard Iteration (or Newton Iteration) can be applied to 


solve the nonlinear system. 


2.4.2  FRONT-FIXING METHOD FOR THE 2D DIFFUSIVE LOGISTIC MODEL WITH RADIAL 


SYMMETRY 


Let us transform the 2D diffusive logistic model with radial symmetry (2.21)-(2.24) 


into a problem with a fixed domain |0, 1]. Under the Landau transformation [12, 34] 


2(t,r) = FO} W(t,z) = U(t,r), (2.27) 


the moving front problem (2.21)-(2.24) turns into 


aw ow D 2G’) 
ee at Tl pg male ee ) 


aw _ 


5, = GW (a— BW), (2.28) 


where t > 0,0 < z <1 and G(t) = H?(t). 
The boundary conditions (2.22) and Stefan condition (2.23) take the forms 


ri 0) =0; W(t,1)=0, t>0, (2.29) 
and 
CO= =e I). aod. (2.30) 


respectively, while the initial conditions (2.24) become: 
G(0)= H3, W(0,z) =W(z) =Up(zHp), O<2z<1. (2.31) 
Conditions (2.24) for the initial function Uo(r) are translated to Wo(z) as follows: 
Wo(z) € C7([0,1]), W5(0)=Wo(1) =0, Wo(z)>0, O<2z<1. (2.32) 


Le 


After the transformation, the new problem is to solve the nonlinear parabolic 
partial differential equations (2.28)-(4.22) in a fixed domain [0,7] x [0,1] for the 
variables (t, z). Let us consider the step size discretization k = At, h = Az =1/M, 
and the mesh points (t", z;), with t? = kn, n > 0, z; = jh,0 <j < M and M is 
the total number of the subintervals of [0,1]. Let us denote the approximate value of 


W(t”, z;) at the mesh point (t”, z;), 


we = Wt", 2;), (2:33) 


and let g” be the approximation of G(t"). Then we consider the forward approxima- 


tion of the time derivatives, 


wrtt_w? aw Gore" 
j Ip fo) 2 wale), 2.34 
and the central approximation of the spatial derivatives, 
wi wha OW.) wh 2ut tw, OW, 
x bees x bo eels 2.35 
oh Oz ( 25); h2 Oz2 ( 25) ( ) 


From (2.34) and (4.25), Equation (2.28) is approximated by: 


Pos we 7 pun — wz + wh. 7 _ 5 grt g” wie — wit 
k h? Z- 2h 
=g"wi(a—bw}), n20, 0<j7<5M—-1. (2.36) 


For the point at 7 = 0, the value w”, is eliminated from the second order discretization 
of the boundary condition (2.29) and (4.22), 


wy aw", 


Transformed Stefan condition (4.20) is discretized using first order forward approxi- 
mation for G’(t) and three points backward spatial approximation of OM it, 1: 


nr 


—* (3why — dw", tw»), n>0. (2.38) 


18 


to preserve accuracy of order O(k) + O(h?). 
Finally, we have 


2Dk ith cess Daas 
he + k(a — bug))wo + apne 


wy"! = (1 


1 < 
wet =arwe yt bfwet+cfui,, n20, 0<j<M-1. 


eet} 
where the coefficients are given by: 


Dk Dk _ 2ypk(4why_y — why») 


n 


a; im gth2 = 2hz;9” Ah2g” ? 
; 2Dk i 

n Dk Dk 25 pk(4wiy_1 — Wi-2) 

Cc; = . 
J grh? 2he.g” 4h?g” 


2.5 LEVEL SET METHOD FOR THE GENERAL 2D MODEL 


The general 2D diffusive logistic model for the density of population of the invasive 


species U(t, x, y) depending on time ¢ and spatial variables (x, y) has the form 


2 2 
M0. OU OU 


op Digger t Gp) =Ua-W), t>0 @y ER), — (2:39) 


together with the boundary condition 
UE, 028) S0,. #0; (2.40) 
and the Stefan condition 
v(t, x,y) = p|VU (2, y)| n(t, 2, y) = —uVU (Ut, 2,y), t>0, (4, y) € ON(t), (2.41) 


where v(t,z,y) and n(t, x,y) are, respectively, the velocity vector of the boundary 
point (x,y) € OQ(t), and the unit outward normal of 2(t) at (x,y) € ON(t). 


The initial conditions are 
92(0) = 2, U(0, 2, y) = Uo(z, y), (x,y) € 2. (2.42) 
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The initial function U(x, y) is assumed to have the following properties: 
Up € C?(29), Up > 0 in M, Uy =0 on OM. (2.43) 


In what follows, to simplify notations, we will use 7(t) to denote the unknown moving 
boundary 0(2(t). Let us recall that the population is distributed in the domain 
9(t), D > 0 is the dispersal rate and the positive parameters a and 0 are the intrinsic 
growth rate and the intranspecific competition rate, respectively. The parameter 
jt > 0 involved in the Stefan condition (2.41) is the proportionality constant between 
the population gradient at the front and the speed of the moving boundary. 
Following the ideas of [18], here we use a level set approach to effectively capture 
the front at each new time step and a finite difference scheme to solve the parabolic 


equation (2.39) everywhere away from the front. 
Step 1. Initialize level set function @(t, x, y). 


We construct a level set function ¢, such that at any time t, the front is equal to 


the zero level set of ¢, i-e., 


T(t) = {(@,y) E20): Olt,a,y) = OF 


Initially, U(0,x,y) = Uo(x,y) and ¢ is set to equal to the signed distance function 


from the front such that @ is negative in (2 and positive in 125, 


Bd; eceelz,, 
p(0, 2, y) _ 0, Le TG; (2.44) 
—-d “Ee Qo. 


where d is the distance from (x,y) to the front. 
The idea behind the level set method is to move ¢ with the correct speed v at 
the front and followed by updating U(t, x,y). The new position of the front is stored 


implicitly in ¢. With this approach, we avoid the difficulties that arise from explicitly 


20 


tracking the front and thus increase the efficiency to deal with complex interfacial 
geometries. 

Given the velocity field v, at which the front moves, we would construct a speed 
function, V(t,z,y), which is a continuous extension of |v(t, x, y)| from the front 7(t) 


over the whole computational domain. The governing equation of ¢ is then given by 
d+ V|V¢| =0. (2.45) 


This equation will move ¢ with the correct speed at the front by assuring that r(t) 
will always coincide with the zero level set of ¢ at time t. 

To see why this is the case, we need to show that with V chosen this way, the set 
T(t) := {(x,y) : (t,x, y) = 0} coincides with 7(¢t) for all t > 0, or equivalently (2.45) 
yields the same equation for the velocity vector vi(t, x,y) at (x,y) € 71(t) whenever 
T(t) coincides with 7(t) at some t > 0 (note that they coincide at time t = 0 by 
assumption). 

Indeed, if 7(t) = 7(t), from @(t, 71(t)) = 0 and (t,x, y) < 0 for (x, y) lying inside 


T(t), we deduce 
d:+Vob-v1 =0, n= V6/|VO| for (x,y) € n(t) = T(E), 
and v, has the same direction as n, the unit outward normal of 7,(t) = T(t), ie., 
vi = V,n for some V; = Vi (t, 2, y) > 0, (x,y) € 1 (t) = T(t). 


These relations yield 


é: +VUi |Vo| = 0 on 7r(t). 


Combining this with (2.45), we obtain 
V, = V for (x,y) € r(t). 


Therefore 


vi = Vn for (a, y) € T(t). 


pal 


By (2.41), we have 


v =V n for (2, y) € T(t), 


and thus we have proved 


vi = Vv for (x,y) € T(t), 
as wanted. 
Step 2. Compute the velocity function V(t, z, y). 


In the algorithm, we compute V(t, z,y) at every grid point. One issue in com- 
puting VU arises from the fact that its approximation is usually in the order O(1) 
at points close to or on the front. By introducing V defined as an extension of 
|v(t,z,y)| = p|VU(t, 2, y)| from (x,y) € T(t), we can avoid unnecessary numerical 
difficulties when we solve Equation (2.45). 

The approximation to VU at 7(t) is based upon approximations to the derivatives 


of U in four coordinate directions to cut down on grid orientation effects. 


woes} ((2) +2) (+) eo 


Each approximation to a derivative of U can be continuously extended away from the 


front by the advection equations 


¢ 


Figure 2.5: Four coordinate directions. 


uz + S(¢bxr)uy = 0, (2.47) 


Ze 


uy + S(doy)uz = 0, (2.48) 
u; + S(bbn)u;, = 0, (2.49) 
ut + S(ddc)ue = 0, (2.50) 


where ut = 0U/Oz, u? = 0U/dy, u® = OU/On and u* = OU/O¢ on r(t). Here J is 
the sign function. 

Equation (2.47) through Equation (2.50) have the effect of continuously extending 
u', u?, u, u* away from the front by advecting these fields in the proper upwind 
direction, and they are used to define V away from r(t). Note that these equations 


will not degrade the value of V on the front because ¢ is zero on T(t), hence, so are 


S(obx), S(Oby), S(bbn) and S(ege). 


Step 3. Update ¢ to be a signed distance function for one time step. 


The computation of the velocity vector at the front is dependent upon the level 
set function ¢. However, by Equation (2.45), the level set function will cease to be 
an exact distance function even after one time step. In order to keep the accuracy of 
n and V, we need to avoid having steep or flat gradients developed in ¢. One way to 
avoid these numerical difficulties is to reinitialize ¢ to be an exact distance function 
from the evolving front 7(t) at each time step. 

Firstly, we update @ by solving (2.45) with V given in Step 2 and initial value 
@ obtained in Step 1. We then follow the same ideas as in [56] to reinitialize ¢. In 
that paper, an algorithm was presented for reinitializing the level set function ¢ to be 
an exact signed distance function from the front. The basic idea behind this method 
is that given a function 9 that is not a distance function, one can evolve it into a 
function @ that is an exact signed distance function from the zero level set of do. This 


can be accomplished by iterating the following equation to a steady state 


ot = 5(¢0)(1 — |V9l), (2.51) 
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where $(0, x,y) = ¢0(z, y) and S again denotes the sign function. As in [56], the sign 


function S is smoothed by the equation 


o 


S-(¢o) = Wee fe 


to avoid numerical difficulties while implemented. 


(2.52) 


By using this approach, we do not need to explicitly find the contour ¢9 = 0, and 
update the values of the front 9 at grid points instead. From Equation (2.51), it is 
clear that the original position of the front will not change, but at points away from 


T(t), @ will be evolved into a distance function. 
Step 4. Update U(t, x, y). 


After moving @ by the correct velocity at the front and then reinitializing @ to 
be very close to an exact signed distance function from r(t) in Step 3, next we 
update U(t, x,y). Updating U(t, x,y) essentially boils down to solving the nonlinear 
parabolic partial difference equation (2.39) over the whole computational domain in 


the following four cases: 


i,j-1 i+f,j-1 


Figure 2.6: Illustration of updating U at four kinds of points. 


e At points away from the front, which means the nearby four grid points are 


all inside the domain {2(t), we solve the nonlinear parabolic partial difference 
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equation (2.39) using a standard five-point stencil. For example, the grid point 


(i+ 1,7) in Figure 2.6. 


e For points near the front T(t), some special care should be taken. We effectively 
capture the front using the level set function @. We can use one-sided different 
sign of @ to incorporate the distances between a point on the front and grid 
points neighboring it in either the vertical or horizontal direction. For example, 
yr € T(t), yp = ((i—1)h, yr) € T(t)we consider two grid points (7,7 +1) and 
(i,j) which border ys. In y-direction, we have y; < yr < yji1. We introduce 


ij 
up — yj = Th =-~—“8__h 
| Pig+ — i,j 

and use U;,;, Uj;-1, Uij;-2, r and U(ys) = 0 to construct extrapolating poly- 
nomial P(y). When we update U;,;, we use a standard five-point stencil by 


employing P(jh) instead of Uj,;41. For the case when front interacts with 2- 


axis, we use the same process in x-direction. 


e For the extreme configuration, where we cannot find enough grid points inside 
the domain to construct interpolating polynomial P, we employ the nearby grid 
points in the domain, and also the intersect points of the front and x and y-axis 


instead, to do Taylor expansion at that grid point to update U. 


e If a grid point lies on the front, we set the value U = 0 at that point (in view 


of (2.40)). For example, we set U;_;,;=0 for the grid point (¢ — 1,7). 


Step 5. Repeat Step 2 through Step 5 to update @ and U for the next time step. 
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2.6 NUMERICAL EXPERIMENTS 


2.6.1 NUMERICAL TESTS FOR FRONT-TRACKING METHOD AND FRONT-FIXING METHOD 


OF 1D MODEL 


Convergence test for front-tracking method of 1D model 

In the logistic diffusion model (2.1)-(2.5) with parameters values (D, 1, a, b, Ho)= 
(0.4,10,1,1,1) and Up = cos(*%"). Here we vary the temporal size with fine spatial 
size. In Table 2.1, the error (both Lz and L,.) and the convergence to the solution of 
front-tracking method is examined, with final time t.,q = 1. The error is computed by 
the difference of the numerical solution with the exact solution. For all the examples 
below when the exact solution is not given, the solution with a fine resolution will be 
considered as reference or "exact" solution. Second-order convergence in the spatial 
dimension can be observed. 


Table 2.1: Accuracy results for front-tracking method of 1D model 


Nx xX Ny LoError | Order | L.Error | Order 
Accuracy test of U of front-tracking method 
61 x 2e06 4.10e-03 4.4e-03 
121 x 2e06 9.40e-04 2.25 9.35e-04 2.14 
241 x 2e06 2.20e-04 2.16 2.10e-04 2.10 
481x2e06 || 4.36e-05 2.36 | 4.08e-05 | 2.33 
961x2e06 || Reference 
Accuracy test of H of front-tracking method 
61 x 2e06 1.68e-01 4.28e-02 
121 x 2e06 2.72e-02 2.63 9.40e-03 2.19 
241 x 2e06 4.4e-03 2.62 2.1e-03 2.15 
481x2e06 || 6.20¢e-04 2.84 | 4.09e-04 | 2.36 
961x2e06 || Reference 


Convergence test of front-fixing method for 1D model 
In the logistic diffusion model (2.1)-(2.5) with parameters values (D, 1, a,b, Ho)= 
(0.4,10,1,1,1) and Up = cos(%). Here we vary the temporal size with fine spatial 
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size. In Table 2.2, the error (both Lz and L,.) and the convergence to the solution of 
front-fixing method is examined, with final time t.,g = 1. As expected, a second-order 
convergence in space can be observed. 


Table 2.2: Accuracy results for front-fixing method of 1D model 


Nz xX Ni LoError | Order | L.Error | Order 
Accuracy test of U of front-fixing method 
21x 1.6384e6 1.40e-03 1.92e-03 
41x 1.6384e6 3.48e-04 1.96 4.83e-04 1.92 
81x 1.6384e6 8.56e-05 2.00 1.19e-04 | 2.02 
161 x1.6384e6 || 2.03e-05 2.07 | 2.84e-05 2.07 
321 x1.6384e6 || Reference 


Accuracy test of H of front-fixing method 


21x 1.6384e6 3.89e-02 7.56e-03 

41x 1.6384e6 7.36e-03 2.17 | 2.00e-03 iiyal 
81x 1.6384e6 1.32e-03 2.40 5.07e-04 1.92 
161 x1.6384e6 2.26¢0-04 2.48 1.22e-04 1.98 
321x1.6384e6 || Reference 


Comparison between front-tracking method and front-fixing method of 1D model 

In Figure 2.7, we use the front-tracking method and front-fixing method to sim- 
ulate the logistic diffusion model (2.1)-(2.5) with parameters values (D, py, a,b, H))= 
(0.4, 10, 1,1,1), Uo = cos() and spatial size h = 0.0125. Figure 2.7 shows that the 
results of front-tracking method and the results of front-fixing method are consistent 


with each other. 


2.6.2 NUMERICAL TESTS FOR FRONT-TRACKING METHOD AND FRONT-FIXING METHOD 


OF 2D MODEL WITH RADIAL SYMMETRY 


Convergence test for front-tracking method of 2D model with radial symmetry 
We consider the 2D logistic diffusion model with radial symmetry (2.21)-(2.24) 
with parameters (D, 1, a,b, Hp)=(0.4, 10,1, 1,0.5) and Up = cos(*}). The system is 


used to test front-tracking method. In Table 2.3, the error (both Lz and L,,) and the 


yas 


—front—tracking for 1d 
‘-« front-fixing for 1d 


u(x,t=1) 


Figure 2.7: A comparison between front-tracking method and front-fixing method of 
1D diffusive logistic model 

order of convergence in space to the solution of front-tracking method is examined, 
with final time 7’ = 0.01. Again second-order convergence in space can be observed. 
Convergence test for front-fixing method of 2D model with radial symmetry 


Table 2.3: Accuracy results for front-tracking method of 2D model with radial sym- 


metry 


N, x N; L,Error | Order | L,.Error | Order 
Accuracy test of U of front-tracking method 
71x 2e04 6.50e-04 2.71¢e-03 
141 x 2e04 1.42e-04 | 2.19 | 5.96e-04 | 2.19 
281 x 2e04 3.24e-05 | 2.14 | 1.85e04 | 2.14 
561 x2e04 6.27e-06 | 2.37 | 2.61le-05 | 2.37 
1121x2e04 || Reference 
Accuracy test of H of front-tracking method 
71x 2e04 3.02¢e-02 5.01e-03 
141 x 2e04 6.75e-03 | 2.16 | 1.07e-03 | 2.23 
281 x 2e04 1.54e-03 | 2.14 | 2.42e-04 | 2.14 
561 x 2e04 3.01e-04 | 2.35 | 4.67e-05 | 2.37 
1121x2e04 || Reference 


We test front-fixing method for solving the 2D logistic diffusion model with ra- 
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dial symmetry (2.21)-(2.24) with parameters (D, 1, a,b, Hp)=(0.4,1,1,1,1) and Up = 
cos(). In Table 2.4, the error (both Lz and L..) and the order of accuracy in 
space of front-fixing method is examined, with final time 7’ = 0.5. As expected, a 
second-order convergence in space can be observed. 


Table 2.4: Accuracy results for front-fixing method of 2D model with radial symmetry 


N, x N, || LeError | Order | L,.Error | Order 
Accuracy test of U of front-fixing method 
21x5e4 || 6.20e-03 8.70e-03 
41x 5e4 1.50e-03 | 2.01 | 2.20e-03 | 2.00 
81x5e4 || 4.00e-04 | 2.07 | 5.00e-04 | 2.07 
161x5e4 || 1.01e-04 2.32 1.00e-04 2.32 
321x5e4 || Reference 
Accuracy test of H of front-fixing method 
21x5e4 1.10e-03 1.90e-03 
41x5e4 || 3.00e-04 | 2.00 5e-05 1.98 
81x 5e4 1.00e-04 | 2.05 le-05 2.05 
161x5e4 || 2.01e-05 2.31 2.02e-06 2.31 
321x5e4 | Reference 


Comparison between front-tracking and front-fixing of 2D radially symmetric model 

In this section, we use front-tracking method and front-fixing method to simulate 
the 2D logistic diffusion model with radial symmetry (2.21)-(2.24) with parameters 
(D, p, a, b, Hp)=(0.4, 10, 1,1,1), Uo = cos(*4}) and spatial size h = 0.00625. Figure 
2.8 reveals that front-tracking method matches well with front-fixing method for the 


2D logistic diffusion model with radial symmetry (2.21)-(2.24). 


2.6.3. NUMERICAL TESTS FOR LEVEL SET METHOD OF 2D MODEL 


Convergence test for level set method of 2D model with radial symmetry 
Here we study the 2D logistic diffusion model (2.39)-(2.43) by using level set 


approach with parameter (D, 14,a,b)=(0.4,10,1,1), 7 is a circle with radius 1 and 


Up = 4cos(n~ ee) 
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si — front-fixing for 2d radial symmetry 
‘»» front-tracking for 2d radial symmetry 


=0.1) 


u(r,t 


Figure 2.8: A comparison between front-tracking method and front-fixing method for 
2D diffusive logistic model with radial symmetry 


For the boundaries of the species, we use the red dotted curve to show the simu- 
lated boundary of the species, the blue circle is introduced to describe to what degree 
the boundary evolves like a circle. The radius of the blue circle is the average distance 
between the intersect points of T(t) with x-axis and y-axis on the boundary and the 
origin, i.e. 


_lVEFP 


Hof (x,y) 

where (x,y) € T(t) are all the intersect points of v(t) with x-axis and y-axis. 

According to [13], the solution of Equation (2.21)-(2.24) is unique and radially 
symmetric. Figure 2.9 shows the evolution of U(t,xz,y) and r(t), where we can see 
that the blue circle matches exactly with the red dotted curve, which means that the 
boundary 7(t) keeps the geometry. And it can be easily observed that U(t,x,y) has 
radial symmetry as Up. 

We focus on the radius of the boundary 7(t), which we denote by H(t). U(t,r) = 
U(t, x,y) is used to learn about the order of accuracy in space of the level set method. 


The convergence test for the solution of u(r) at T = 0.1 and the front H(t) can 
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Figure 2.9: Numerical simulation of the density of population u(t, x,y) and moving 
boundary 7(t) with initial domain (29 a disk in 2D using level set method. 


be observed from Figure 2.10 with different space sizes h = 0.025,h = 0.0125,h = 


0.00625, h = 0.003125, and the results are compared to the results of front tracking 


method with the same initial setup and step size h = 0.003125. Figure 2.10 shows 


that the results of level set method agree very nicely with the results of front tracking 


method, which means the three methods are consistent with each other. 


In Table 2.5, the error (both Ly and L.,) and the order of convergence to the 


solution of level set method is examined, with final time T = 0.1. It reveals that the 


convergence orders for both the solution u and the front H(t) are between 1 and 2. 


2.6.4 NUMERICAL DICHOTOMY: SPREADING VERSUS VANISHING 


Example 1. In the diffusion logistic model with a free boundary (2.1)-(2.5), we take 


the parameters values (D, p, a,b, Hp)=(0.5,5,2,1,1) and Up = cos(3,-). Figure 2.11 


shows the spreading behavior under the condition Hp = 1 > L = 0.785. 


set h=0.025 
set h=0.0125, 


st h=0.00625, 


I 

le 
le 
{ t h=0.003125 


=: level set 
fronttracking h=0.003125 a 
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t 


Figure 2.10: Convergence test for level set method of 2D model with radial symmetry. 


Table 2.5: Accuracy results for level set method of 2D model with radial symmetry 


Nz X Ny X Nz LoError | Order | L.Error | Order 
Accuracy test of U of level set method 
29x 29x 160 5.58e-03 9.29e-03 
57x57x640 3.06e-03 | 0.86 | 5.01le-03 | 0.89 
113x113 x 2560 1.40e-03 1.13 2.26e-03 1.15 
225 22510240 || 4.84¢e-04 1.54 7.79e-04 1.54 
449 x449x 40960 | Reference 
Accuracy test of H of level set method 
29x 29x 160 4.19e-02 5.84e-02 
57x57x640 2.01e-02 1.06 | 2.76e-02 | 1.08 
113x113 x 2560 8.70e-03 1.20 1.19e-02 1.22 
225 X 22510240 || 2.91e-03 1.57 3.92e-03 1.60 
449 x449x 40960 | Reference 


Example 2. In the diffusion logistic model with a free boundary (2.1)-(2.5) with 


the parameters values (D, py, a,b, Ho)=(1,5,1,1,0.496) and Up = cos(#*). 


Sho Figure 


2.12 shows the spreading behavior occurs even if Hp = 0.496 < L = 1.571. 
Example 3. In the diffusion logistic model with a free boundary (2.1)-(2.5), 
parameters values are set as (D, py, a,b, Ho)=(1,5,1,1,0.496) and Up = sCOsla 


As in this example, we keep all the parameters values same with Example 2 except 


that we decrease the initial value Up, then the vanishing behavior occurs under the 


32 


The density of population at t=5 


The free boundary 


Figure 2.11: Numerical simulation of spreading under condition Ho > L. 


The density of the population at t=20 


The free boundary 


Figure 2.12: Numerical simulation of spreading under condition Ho < L. 


condition Hy = 0.496 < L = 1.571. Figure 2.13 shows the vanishing behavior under 
the condition Hp = 0.496 < L = 1.571. 


9 x10" The density of the population at t=20 The free boundary 
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Figure 2.13: Numerical simulation of vanishing under condition Ho < L. 
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2.6.5 NUMERICAL EXPERIMENTS FOR LEVEL SET METHOD OF 2D MODEL WITH 


DIFFERENT INITIAL CONFIGURATION 


Example 4. In the 2D logistic diffusion model (2.39)-(2.43) with parameters 
(D, uw, a,b) = (4,10,1,1), the initial boundary 7 is set to be an equilateral triangle 
which centers at the origin point (0,0) with side-length 1. The initial value uo(z, y) 


and the initial level set function ¢o(x, y) are set as following 


Uo(x, y) = 
0 (x,y) € 96. 
(V3a—-y+ +e) (-V3a—-y+ ) 
—min(S = A r Y, 5) vs ’ 2 vs ’ (x, y) 2 2o, 
do(x, y) — 0 (x,y) = TO, 
: |\V3a-y+| |-V3e-yt++| Ac 
min(|¥3 = A t yl, 5, 5), (at, y) € 2. 


For the boundaries of the species, we use the red dotted curve to show the sim- 
ulated boundary of the species, the green triangle represents the initial boundary. 
Figure 2.14 shows the simulation of the evolvement of the species and moving bound- 
aries along time with an equilateral triangle as the initial boundary. From Figure 
2.14, we can see that the red dotted curve evolves into a circle, and then propagate 


as a circle, which also agrees with the theoretical results [16]. 


Example 5. In the 2D logistic diffusion model (2.39)-(2.43) with parameters 
(D, uw, a,b) = (5,10,1,1), the initial boundary of the species 7) is a rectangle with 
length=0.6 and width=0.5, centered at (0,0). And the initial function uo(z,y) and 


the initial level set function ¢0(z, y) are set as following 


200(0.5 — z)(0.5 +.2)(0.6—y)(0.6+y), (2,y) € Mp, 
ULL, Y) = 
0 (x,y) € 926. 
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Figure 2.14: Numerical simulation of the density of population u(t, x,y) and moving 
boundary 7(t) with initial domain an equilateral triangle in 2D using level set 
method. 
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—min(0.5 — |x|,0.6 — |y]), (a, y) € 2, 
oo(2,y) = 4 0 (x,y) € To, 
min(|0.5 — |a||,|0.6 — |yll), (wy) € 2%. 

For the boundaries of the species, we use the red dotted curve to show the sim- 
ulated boundary of the species, the green rectangle represents the initial boundary. 
Figure 2.15 shows the spreading of u and moving boundary along time with a rectan- 
gle as the initial boundary. It indicates that the boundary evolves into a circle, and 
then propagates like a circle. 

Example 6. Here we test the level set method for solving (2.39)-(2.43) with two 
other different initial domain setup: annulus (Figure 2.16) and annulus with a cut 
(Figure 2.17). For the boundaries of the species, we use the red dotted curve to show 
the simulated boundary of the species, the green dotted curve represents the initial 
boundary. For all two different cases, the front will asymptotically evolve into circles 


that correlates exactly with theoretical results [16]. 


2.6.6 NUMERICAL TEST FOR LEVEL SET METHOD OF 2D 


ADVECTION-REACTION- DIFFUSION MODEL 


We consider a 2D Advection-Reaction-Diffusion (ARD) model with free boundary of 


the form 

OU OU PU OU ou 

——D) | = a . 
a Daa + goa) t UG, + G,) = Ula BW), t>0, (eu) € M(H), (253) 


together with the boundary condition 
Ul O)).=0. 4 >0, (2.54) 
the Stefan condition 


v(t,v,y) = wIVU(E, ©, y)|n(t,2,y) = VUE, @,y), t>0, (,y) € AQ), (2.55) 
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Figure 2.15: Numerical simulation of the density of population u(t, x,y) and moving 
boundary 7(t) with initial domain {2 a rectangle in 2D using level set method. 
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Figure 2.16: Numerical simulation of the density of population u(t, x,y) and moving 


boundary 7(t) with initial domain 2 an annulus in 2D using level set method. 
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Figure 2.17: Numerical simulation of the density of population u(t, x,y) and moving 
boundary 7(t) with initial domain 9 an annulus with a cut in 2D using level set 
method. 


where v(t,z,y) and n(t, x,y) are, respectively, the velocity vector of the boundary 


point (x,y) € O(t), and the unit outward normal of ((t) at (x,y) € OM(t). 


38 


The initial conditions are 
92(0) = 2, U(0, x,y) — Uo(z, y), (x,y) e 2. (2.56) 


The initial function Up(x,y) is assumed to satisfy (2.43). 

In (2.53), the advection term 3(94 + o) is in the north-east direction. We may 
think of (2.53) as describing the spreading of a flying insect species U affected by 
wind blowing to the north-east direction during the spreading process. 

In the 2D ARD model (2.53)-(2.56) with parameters (D, 1, a,b, 3) = (10, 10, 1, 1,50), 
the initial boundary of the species 79 is a circle with radius equals 1.5, centered at 


(0,0). And the initial value uo(zx, y) and the initial level set function ¢o(x, y) are set 


as follows 


6cos( V oa ym), (2, y) € 2, 
0, (x,y) € 2%; 


do(x,y) = —(0.5 — \/a? + y?). 


Figure 2.18 shows the spreading of U and the moving boundary of this ARD model 


uo(2, y) a 


as time increases, where the red dotted curve represents the simulated boundary of 
the species. In order to clearly reveal the effect of the advection in the model, the 
initial boundary is indicated in the graph by the green circle; the free boundary clearly 
expands faster in the north-east direction and slower in the south-west direction, due 


to the advection in the north-east direction. 
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Figure 2.18: Numerical simulation of u(t, xz, y) and r(t) for initial boundary a circle 
in 2D with an advection term using level set method. 
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CHAPTER 3 
NUMERICAL METHODS FOR A DIFFUSIVE COMPETITION 


SYSTEM WITH FREE BOUNDARIES! 


''S. Liu and X. Liu. Mathematics, 6 (2018), 72. 
Reprinted here with permission of publisher, 11/13/2019. 
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3.1 ABSTRACT 


Numerical methods of the spreading behaviors of two invasive species modeled by 
a Lotka-Volterra diffusive competition system with two free boundaries, which are 
natural extensions of the corresponding free boundary problems of reaction-diffusion 
system, are studies in this chapter. We extent front-tracking method, front-fixing 
method and level set method to the two-species competition system with free bound- 


aries. Results will show these methods are extremely successful. 


3.2. INTRODUCTION 


In [17], Dr. Du and Dr. Wu investigate the spreading behavior of two invasive species 


modeled by a Lotka-Volterra diffusive competition system with two free boundaries 


in RX (N > 1) with spherical symmetry. This system models the invasion of an 


empty favorable habitat, by two competing species U and V, each obeying a logistic 
growth equation and invading the environment through their own free boundaries: 
the species U has a spreading front at r = 5S;(t), while the species V’s spreading 
front is at r = S)(t). The competition-diffusion model for the density of population 
of the competing species U(t,r) and V(t,r) depending on time t and spatial location 


r states as follows: 


U,—D,AU =y,U1-U-kV), t>0, 0<r< S,(t), 
V,-— DAV =yV1-V—-KWU), t>0, O<r< S(t), 
U,(t,0) = V,(¢,0) =0, ¢20, 

(P) 4 Si(t) = —mU,(t, Sit), S3(t) = —HeV-(t, Sa(t)),  t 2 0, 
VERS). or ee Sie. VEO =O for re Sai aad, 
U(0,r) =Uo(r), V(0,r)=Vo(r), for r € [0, 00), 

SOS oe Se. S50) S55 0: 


They show that, for the weak-strong competition case (0 < k < 1 < h), under 
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suitable assumptions, both species in the system can successfully spread into the 
available environment, but their spreading speeds are different, and their population 
masses tend to segregate, with the slower spreading competitor having its population 
concentrating on an expanding ball, say B,;, and the faster spreading competitor 
concentrating on a spherical shell outside B; that disappears to infinity as time goes 
to infinity. 

We aim to treat two invasive species modeled by a Lotka-Volterra diffusive com- 
petition system with two free boundaries in one dimension and higher dimensions 
numerically, which can be a continuation and numerical complement of [17, 60]. We 
propose front-tracking method and front-fixing method in one dimension, we also 
introduce level set method for general two-dimensional case. Illustration with numer- 
ical examples of the spreading behavior of numerical solutions of two invasive species 
modeled by a Lotka-Volterra diffusive competition system with two free boundaries 


are given. 


3.3 NUMERICAL METHODS FOR 1D TWO-SPECIES COMPETITION-DIFFUSION 


MODEL 


Here we study the spreading behavior of two competing species in 1 dimension de- 


scribed by the following free boundary problem: 


U,— Diy, = WU -U-kiV), t>0, O<2< S(t), (3.1) 
Vi — DoVex = V2V(1-V— KU), t>0, O0<2< S(t), (3.2) 
CE = VO SOs 20! (3.3) 
Si(t) = —mU z(t, Sit), 93(t) = —HeVe(t, Sa(t)), t 2 0, (3.4) 


Cia s0: pore Sit). Vise 0° fore 05 (0). ace Oy ~ 13:5) 
U(0,2) =Uo(x), V(0,2) =Vo(x), for x € (0,00), (3.6) 


SO) SS) 0. 85/0) 385 0. (3.7) 
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where U(t, x) and V(t, x) represent the problem densities of the two competing species 
at location x and time t. 7, and 72 are the intrinsic rate of increase of species U and 
V separately. AK , represents the competition effect of species V on species U, and 
Ky represents the competition effect of species U on species V. We are interested in 


investigating the long-term behavior of the free boundary problem numerically. 


3.3.1 FRONT-TRACKING METHOD FOR 1D TWO-SPECIES COMPETITION-DIFFUSION 


MODEL 


The problem lies in solving the nonlinear parabolic partial differential equations (3.1)- 
(3.7) in the fixed computational domain [0,7] x [0, Z] for the variables (t, 2). Let us 
consider the step size discretization k = At, h = Ax = L/M, and the mesh points 
(t”, x), witht” = kn,n > 0,2; = jh,0 <j < M and M is the number of subintervals 
in [0, ZL]. Let us denote the approximate value of U(t", x;) and the approximate value 


of V(t", z;) at the mesh point (t”, x;) 


u; ~U(E", 23), vp & V(t", 25). 


Step 1. Track the position of the moving front 5;(t) of Species U. 


According to the Stefan condition 
Si (t) = —p4U Z(t, Si(t)), t > 0, 


we consider using the central approximation of the spatial derivatives to approximate 


&'(t, $;(t)), which can be divided into the following four cases: 


SP -2j 
h 


1. When 2; < S? < 241, 1 = 2,3,...,.M4 —1, denoting d= , we consider the 
symmetric point of z;_; respect to the position Sj’, which is denoted by z;_1. Espe- 
cially when S? = 2, %j-1 = X41. We use the Lagrange extrapolation to construct 


polynomial P” from the value of d, h, u?_,, u?_,, u? and S$”, thus at %;_1, we use the 


value of P¥ at %;_ instead of u(t”, %;_1), 


44 


OU = P?(%;_1) — Un 4 


eae Gaol =2,.3,...,.M-—-1 

Aa | ) 1) 2(1+ d)h ) a ae) ) 

Zo Lig? Liet Xi Sy Tig th 
TH THT” 


(l+dn (1+a@h 
Figure 3.1: Illustration of how to evaluate Mt, Si(t)) at t” when x; < SP? < x41, 


(=2,3..M—1. 


2. When 0 = x < Sj) < x, the central scheme approximation of the spatial 


derivatives to approximate St, Si(t)) involves the fictitious value u”, at the point 
(t", —h). The value u”, can be estimated from the second-order discretization of the 


boundary condition (3.3), 
Cy 
aT =0 (3.8) 
which implies that wu”, = uf = 0. It is obvious that all the values of u? on the 
grid points are equal to 0 except uj. Numerically, we take S}(t) = 0, and it can 
be explained that the species is only located inside one grid mesh. The simulation 
should stop here indicating that a more refined mesh is needed. 


Stra 


3. When 2; < S} < x, denoting d= 2 7 Let us first consider the symmetric 
point of x9 respect to the position S?, which is denoted by Zp. Then we consider the 
value of u”, = u%, and use the Lagrange interpolation to construct polynomial P” 


from the value of h, d, u”,, u?, uj’ and S?. Then at Zo, we use the value of P” at Zo 


instead of u(t”, Zo). 


a ree XO Ey SP Xo Ci een <M 
0. 
(1+d)h (1+d)h 


Figure 3.2: Illustration of how to evaluate 24(t, S;(t)) at t” when x1 < S? < 2. 


4. When S} = xy, it implies that the spreading of the populations already goes 


out of the computational domain [0, Z], and the simulation should stop here. 
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Step 2. Repeat Step1 to S(t) to track the position of S(t) of species V. 


Step 3. Update the value of U(t"'', x;) and V(t""!, 2). 


1. When x; = S?*" and Le Sot! we know that the moving front of species 
U and the moving front of species V are located on grid points. Setting u?*' = 0, 
uy? = 0, for? =i+1,1+2,...,M and v?** =0, ut! = 0, form =j+1,j+2,...,M. 
We consider the central approximation of the spatial derivatives U,, at x), for 1 = 
0,1,2,...,2-—1, and the central approximation of the spatial derivatives V;. at 2m, for 


m =0,1,2,...,7 — 1, where U and V are updated by the backward Euler in time 


ntl n n+l n+l n+l 
Up — Uy Upp — 2up + uy n+1 n+1 n+1 
a ane =D, ip tyup (1 — uy — Kyu”). 
+1 n n+1 n+l n+l 
ue’ — y Vv —2v +v 
m m _ m-1 m m+1 n+1 n+1 n+1 
eat asic Dy 12 t yur (1l—uh — Keun). 


where / = 0,1,..4 1 and m = 0,1,...j — 1. Then we use the Picard Iteration (or 


Newton Iteration) to solve the nonlinear system. 


Z gvti_ : 
2. When xz; < Ser <i ONG! Shae < %j41, denoting Ry = +5 “ and 
m+1_ . é 
he= 82 ; *: we use the Lagrange extrapolation to construct polynomial P/ from 


the value of h, Ry, us, ut}, ult and Si'* and polynomial P} from the value of 


h, Ra, unt unt, ee and St. Then at x41 and x;41, we use the value of P/ at 


Ti41 instead of uw’! and the value of P/ at x41 instead of ee For the solution u at 


x, for 1 = 0,1,2,...,2 — 1, a standard central approximation in space with backward 


Euler in time will be employed. u/t! = 0, for] =i+1,...M. For the solution v at 
x, for 1 = 0,1,2,...,7 — 1, a standard central approximation in space with backward 


Euler in time will be employed. v;'*t = 0, for] = 7 +1,...M. U and V are updated 
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by the backward Euler in time 


ntl n n+l n+l n+l 

te a Uy — au 

a = dD, 72 oe yup (1 =qpr = Kyu"). 
n+l n n+1 n+1 L 

Uw — Us Ue, — Qu” + Py (ein 
t ; ae dD, a a = 1 ( ) yur (1 _ yett = Kyu?tt?). 
n+l n n+1 n+1 n+1 

Um Um _ Um-=1 2Um T Um+i1 nt+1 n+1 n+1 

—— Dy 12 youn (L— up’ — Kour). 
n+l _ ayn n+1 n+1 LD fas, 

U5 U; = pe 20; ay ee (xj41) yu (1 —yntl _ Kyu"*?) 

k h2 J J J 


where |! = 0,1,... —1 and m= 0,1,...j — 1. Picard Iteration (or Newton Iteration) 
will be applied to solve the nonlinear system. 


3. When x; = S?** and eS Spite £741, then we know that U(t"t!, x;) = 0. Let 


ult! = 0, ult! =0, for] =i+1,i+2,...,M. We consider the central approximation 
n+1_ a. 
of the spatial derivatives U,, at 2;, for 1 = 0,1,...,i —1. Denoting Ra = 53 A) We 


use the Lagrange extrapolation to construct polynomial P/ from the value of h, Ro, 


unty, ort, ort? and S3*!. Then at 241, we use the value of Py at 2;41 instead of 


DAN where U and V is updated by the backward Euler in time 


n+l n n+1 n+1 n+1 
Uy —U Uy, — 2u;,"° +u 
i ee Cee Di l—1 l [+1 ova = yet K ort) 
k h? 
n+l ayn n+l n+l n+l 
Um Um men Um=1 2Um + Um-+i \ cee el > ttl K ee) 
; = 12 h2 r Y2Um Um QU, )- 
n+l n n+1 n+l1 iL 
Uv, — Ur vet, — Que") + Py (2; 
J i Spy, tot 2 (Rie) yt — yt Kaul) 
k - fe TP285 j a ae 


where / = 0,1,..4-—1 and m = 0,1,...j — 1. Then we use the Picard Iteration (or 


Newton Iteration) to solve the nonlinear system. 


m+1 


f Ss a 
4. When a; < St*’ < a4, and a etl denoting Ry = 2 7 “we use the 


Lagrange extrapolation to construct polynomial P/ from the value of h, Ry, u%*y, 


ubt! ut! Then at x41, we use the value of P/ at z;4, instead of ult!. For the 
solution wu at x), for 1 = 0,1,...,2 —1, a standard central approximation in space with 
backward Euler in time will be employed. uj't? = 0, for] =i+1,i+2,...,M. For the 
solution v at x;, for 1 = 0,1,...,7 — 1, a standard central approximation in space with 


backward Euler in time will be employed. v;’t' = 0, for 1 = 7 +1,j +2,...,M, where 


AT 


U and V are updated by the backward Euler in time 


I n n+1 n+1 n+1 
U) — Uy Ui-1 2uj ae Ul41 n+1 n+ n+1 
E = dD, he 1U) (1 —U Kyu, ) 
n+1 n n+1 n+1 LD 
Uz Ui; wey — 2ug" + Pi‘(tis1) n+1 n+1 n+1 
i =D, 72 byyure (1 — ui Ky") 
n+l _ yn n+l n+l. ,,n+1 
Um Um ns Um—1 2Uim Um-+1 Betty _ gyntl _ K ian) 
ke = L’2 he Y2Um Um Quin): 


where | = 0,1,..4 -—1 and m= 0,1,...j — 1. Picard Iteration (or Newton Iteration) 


will be applied to solve the nonlinear system. 


3.3.2. FRONT-FIXING METHOD FOR 1D TWO-SPECIES COMPETITION-DIFFUSION MODEL 


Here we consider transforming Equation (3.1) and Equation (3.2) into problems with 


a fixed domain |0, 1] separately. 
Step 1. Update the front of S(t) and the value of U by front fixing method. 


Let us transform Equation (3.1) into a problem with a fixed domain [0, 1] using 


the Landau transformation |12, 34] 


Ce ea M(t,y) =U(t,x), WtoM(t,y) = V(t,2). 


x 
Si(t)’ 
Then Equation (3.1) turns into the form: 


OM yOM O*>M 
2 Oy : Oy? 


= H(t)y,M(1— M — K,WtoM), (3.9) 


where ¢ > 0,0 < y <1 and H(t) = S?(t). 
Boundary conditions (3.3) take the forms: 
OM 
o“ (t,0)=0, M(t,1)=0, t>0, 
Oy 
and Stefan condition (3.4) is transformed into 
OM 
H'(t) = —2u,——(t,1), t>0, (3.10) 
Oy 
respectively, while the initial conditions (3.7) become: 


H(0) = (S?)?, M(0,y) = Mo(y) = Uo(yS?), O<y<1. (3.11) 
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Conditions (3.6) for the initial function Up(x) are translated to Mo(y) as follows: 
Mo(y) € C7((0, 1), Mo(0) = Mo(1) =0, Moly) >0, OSy<l. (3.12) 


After the transformation, the new problem has been changed to solve the nonlinear 
parabolic partial differential equations (3.9) in the fixed domain [0,7] x [0,1] for the 
variables (t,y). Let us consider the step size discretization k = At, h = Ay = 1/L, 
and the mesh points (t”, y;), with ’ = kn, n > 0, y; = jh, 0 <j < Land L is the 
number of subintervals of [0,1]. We denote the approximate value of M(t”, y;) at the 


mesh point (t”,y;), 
m” = M(t",y;), wtoM? » WtoM(t", y;). 


and let H” be the approximation of H(t”). Considering the forward approximation 


of the time derivatives, 


1s Jj 


ntl _m™ OM ah etl yn 
ki ~ ap aD ki 


~ H(t"), (3.13) 


and the central approximation of the spatial derivatives, 


n fos n n =. n n 2 
Mey, — m1, OM my 1 — 2p +m, OOM 


i ~ is ; z ~ " i). . 
ah a Oy (t ah) h2 Oy? (t £Y5) (3 14) 


Equation (3.9) is then approximated by 


met! _ mn 


nod J 
ss k 2 2h k 


em, —m” n+1 _ FTn tty ae ie mh 
Yj Miz — M51 A my — 2mz + Mi, 


or 72 


(3.15) 
= H"ymj(1—m} —wtoM?), n20,0<j<M-1. 


As usual, we assume that the Equation (3.15) can also be approximated at j = 0. 
Equation (3.15) written for 7 = 0 involves the fictitious value m”, at the point 
(t",—h). The value m", is eliminated from the discretization of the boundary and 


initial condition (3.11) and (3.12), 


Transformed Stefan condition (3.10) is discretized using first order forward approxi- 


mation for H’(t) and three points backward spatial approximation of a 1): 


Hike 


Att = H™ = 
h 


(8m5, —4mi,_1 + miy_2), 220 


to preserve accuracy of order O(k) + O(h?). 
Finally, we have 


2D,\k Dyk ‘a 
The +ky(1 — mg — Kiwe))ms + 2am 


m5" = asmy_, + imi + mi, n20, O<jSM—1. 


mst = (1 - 


mS; 
where the coefficients are given by 


n Dk yjlak(4miy_4 —m'y_9) 


Oy Fe Ah2H” , 
2Dik 

b= ere aa + ky (1 — mi — KywtoM;), 

n Dik yjbaik(Amiy_ 1 — Miz») 

5 Hp? 4h? H" 


Step 2. Update the front S(t) and the value of V by front fixing method. 


Let us transform Equation (3.2) into a problem with a fixed domain [0, 1] using 


the Landau transformation [12],[{34] 


A a a Wlh2) = Vite). MioW (iz) =U ta): 


a 
So(t)’ 
Then the Equation (3.2) turns into the form: 


ow ., zoW _ &w 
GOS - O55, — eae = GOnRW(L- W - K2Mtow), (3.16) 
where t > 0, 0 < z < 1 and G(t) = S2(t). 


Boundary conditions (3.3) take the forms 


OW 


sz — (t, 0) = 0 W(t,1)=0, t>0 
Dz (9) ’ (g,1) ’ > UV, 
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and Stefan condition (3.4) turns into 


Gis oT (t 1), 0; (3207) 


respectively, while the initial conditions (3.7) become 
G(0) = (S)’, W(0,z) = Wo(z) =VolzS2), O<2<1. (3.18) 
Conditions (3.6) for the initial function Up(x) are translated to Wo(z) as follows: 
Wo(z) € C?([0,1]), W5(0) =Wo(1) =0, Wil(z)>0, O<2z<1. (3.19) 


After the transformation, the new problem lies in solving the nonlinear parabolic 
partial differential equations (3.16) in the fixed domain [0, 7] x [0, 1] for the variables 
(t,z). Let us consider the step size discretization k = At, h = Az = 1/L, and the 
mesh points (t", z;), with ’ =kn, n >0, 2; = jh, 0<j < Land L is the number 
of subinterval of [0,1]. We denote the approximate value of W(t", z;) at the mesh 
point (2" 27); 

w) © Wt", 2;),  mtoW; ~ MtoW(t", z;). 


and let G” be the approximation of G(t"). Let us consider the forward approximation 


of the time derivatives, 


oes WF an (t”, x), a AG Gy, (3.20) 
and the central approximation of the spatial derivatives, 
mia rw an Gz) i a mae ot (es (3.21) 
From (3.20) and (3.21), Equation (3.16) is approximated by: 
nee = uF ay Wh — Wha Get - Sl A fea LE 


k 2 2h k he 


= G"youwi (1 =e KymtoW;'), n>0,0<j<M-1. 
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As usual, we again assume that the Equation (3.22) can be also approximated at 
j = 0. Equation (3.22) written for 7 = 0 involves the fictitious value w”, at the point 
(t",—h). The value w", is eliminated from the discretization of the boundary and 
initial condition (3.18) and (3.19), 


wy — we 3 
ae) eames wr, =0, n>0. 
Transformed Stefan condition (3.17) is discretized using first order forward approxi- 
mation for G’(t) and three points backward spatial approximation of 2 wit, 1): 


Hak 


Grrl = G” = 
h 


to preserve accuracy of order O(k) + O(h?). 


Finally, we have 


2Dok , Dok 
Gny2 + ky2(1 — wo — KomtoW;'))we + 2 Gn wl 


went = Ajwj_,+ Beuy t+ Ceuy, 220, 0<7<5M-1. 


wort = (i= wi 


where the coefficients are given by 


a. Dok 7 Z;pok(4wiy_1 — wiy_a) 
J Grh2 Ah2GQr ’ 


By =1- Gnpe + ky2(1 — w? — KgwtoM?), 
Cr= Dok | ZjMok(4wiy_1 =: wu-2) 
j Gnrp2 - Ah2GQr : 


Step 3. Update WtoM(t", y;) with the front information G’*! and H”*1. 


1. When y; = 1/ _ ;, we know that WtoM(t"*!, y;) = 0. Setting wtoM?'t! = 0, 


wtoM?*'= 0, for! =i+1,i4+2,...,M, we consider the central approximation of the 


spatial derivatives & Oye at y;, for 7 = 0,1,...,i-1, where WtoM is updated by the 


backward Euler in time 


2Dik Dik 
rae + ky(1 — wtoM? — Kim¢))wtoM> + ant 


wtoM?** = = a-wtoM;, + bjwtoM; + c;wtoM7,, n20, O<j<i-l. 


wtoM"' = (1- 


wtoM?’. 
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n+ . nm+1 n+1_ A 
2. When y; < ot < Yj411, denoting R = 5 es we use the Lagrange 


extrapolation to construct polynomial P” from the value of h, R, wtoM?>', wtoM?'4', 


wtoM}"*! and aa then we use the value of P¥ at yj. instead of wtoMfH!. 


2D k 7 Dik 7 
aaa + ky(1 — wtoMy — Kimg))wtoMg + Daan wtoM?. 


wtoM?*? = ai jwtoMe ot bewtoM? + cj ™utoM” jie 20 Ua gpet— 1. 


wtoMf*! = (1- 


wtoM?*" = a?wtoM? , + bPwtoM? + cP P§(y441). 


3. When 4/ oo > 1, we consider the central approximation of the spatial derivatives 


a Oye at y;, for 7 = 0,1,...,M —1. When updating WtoM(t"*", ya), we know 


that W(t", ,/ £5) = WtoM (t"*}, yz), we approximate W(t"t?, ,/ #52) instead of 


approximating WtoM (t"*!, yay). Suppose z; < ort < 241, we use the Lagrange 


extrapolation to construct polynomial P from the value of 2;, 2-1, 2:41, W/",, W,”, 


and W,, thus WtoM (t"*1, yw) = W(t", i) P(g 


Gr 


Yo yy peeved he. Oe Yu ixl 
0) cAI Zl 
) eG st 


Figure 3.3: Updating WtoM (t"*', yy,) in U-system by evaluating W(t"™?, ) in 


V-system instead when case eae 


Step 4. Update MtoW(t”, z;) with the front information G’*! and H”*?. 


1. When z = fen, we know that MtoW(t"*", z;) = 0. Setting mtoW?*? = 
0, mtoW;"*! = 0, for 1 =i+1,i+ 2,...,.M, we consider the central approximation of 


e Mew at z;, for j =0,1,...,i—1, where M is updated by the 


the spatial derivatives 
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backward Euler 


2D2k Dok 
Gnj2 + ky2(1 — mtoWs — KowG))mtoW 7 + 2 ange aya 


mtow?*? = AimtoW;", + B?mtoW; + CimtoW;,, n20, O<jS<i-l. 


mtowst! = (1— mtoW;,’. 


4 H+1/Gnrtl_g, 
2. When 2; < \/ Gan < 241, denoting R = ~—_{~—_—_, we use the Lagrange ex- 


trapolation to construct polynomial P” from the value of h, R, mtoW?S', mtoW;"4', 


mtoW;"*" and fee . Using the value of P” at 2; instead of mtoW ss", then we 
get 


2D2k Dok 
aa + kyo(1 — mtoW — KowG))mtoWs + on nf 


mtow?*? = AimtoW;", + B?mtoW; + C7mtoW;,, n20, O<jS<i-l. 


mtoWgt! = (1- mtoW}. 


mtoW?t! = A®mtoW? , + BPmtoW? + CPP" (z:41). 


3. When 4/ oe > 1, we consider the central approximation of the spatial derivatives 


& oMeOW at z;, for j = 0,1,2,...,M —1. When updating MtoW(t""", zy), we know 


that M(t"*!, ,/ S95) = MtoW(t"*!, 2,1), let’s approximate M(¢t"*+, \/ $5) instead 


er, 2 


of approximating MtoW ( mu). Suppose y; < ot < yi41, we use the Lagrange 


extrapolation to construct polynomial P from the value of y;, yi-1, Yis1, M71, M7, 


and Mf, then MtoW(t"*?, zy) = M(t"*1, \/ S55) » P(,/ S55). 


Figure 3.4: Updating MEE zy) in V-system by evaluating M (t+, / G5) in 
as di. 


Hrt+ 


eat 


U-system instead when 
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3.4 LEVEL SET METHOD FOR 2D TWO-SPECIES COMPETITION-DIFFUSION 


MODEL 


A general 2D two-species competition-diffusive model for the densities of population 
of the species U(t, x, y) and V(t, x, y) depending on time ¢ and spatial variable (x, y) 


states as follows: 


OU CU. OF 
a 7 Pilg alt, Op =yU(1-U—k,V), t> 0, (x,y) € i(¢), (3.23) 
OV Ove ev 
nO 7) epee alle 
at (a2 Oy? 


J=yV(1-V-KU), t>0, (2,y)€ 2o(t), (3.24) 
together with the boundary conditions 

UG nit) y = 0, Viz =, “te 0), (3.25) 
and the Stefan conditions 


vilt,2,y) = w1|VU(t, 2, y)| ni (t, 2, y) = —uVU(t,2,y), t>0, (3.26) 

vo(t, 2, y) = pe|VV(t, 2, y)| na(t, 2, y) = —uwVV(t,z,y), t>0, (3.27) 

where vi(t,z,y) and n,(t, x,y) are, respectively, the velocity vector of the boundary 

point (x,y) € O2,(t), and the unit outward normal of (2,(t) at (x,y) € O22, (t), 

vo(t, x,y) and na(t, x,y) are, respectively, the velocity vector of the boundary point 
(x,y) € O922(t), and the unit outward normal of 22(t) at (xz, y) € O29(t). 

The initial conditions are 
U(0, x,y) = Uo(z,y), (x,y) = 9,(0), (3.28) 


V0, 259) = Vo(a, y), (x,y) S §23(0). (3.29) 
The initial functions Up(x,y) and Vo(z, y) satisfy the following properties: 


Uo(z,y) € C*(21(0)), Ug(0) = Uo(7i(0)) = 0, Uolz,y) > 0, (wy) € 2(0). 


Vo(x, y) € C?(22(0)), V5 (0) = Vo(72(0)) = 0, Vo(x, y) > 0, (e, y) € §22(0). 
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Here 7,(t) and 79(t) are the unknown moving boundaries of two species U(t, x,y) and 
V(t,x,y) such that the population are distributed in the domain 2\(t) and 23(t) 
separately. D; > 0 and Dz > 0 are the dispersal rates. The parameters jz; > 0 and 
[lg > 0 involved in the Stefan conditions (3.26) and (3.27) are the proportionality 
constant between the population gradient at the front and the speed of the moving 
boundary of two species U(t, x,y) and V(t, x,y) respectively. Following the ideas of 
[18], here we use a level set approach to effectively capture the front at each new time 
step and a finite difference scheme to solve the parabolic equation everywhere away 


from the front. 
Step 1. Initialize level set function ¢1(t, x,y) for species U. 


We construct a level set function ¢,, such that at any time t, the front 7,(t) is 


equal to the zero level set of 1, i-e., 


Ti(t) = {(a,y) € A(t) : dilt, x,y) = OF 


Initially, U(0, x,y) = Uo(x, y) and ¢, is set equal to the signed distance function from 


the front of species U such that ¢, is negative in (2;(0) and positive in 2¢(0), 


+d, x € N§(0), 
b1(0,2,y) = 0, x €7,(0), 
—d x €92,(0). 
where d is the distance from (x,y) to the front 7(t). 

The idea behind the level set method is to move ¢; with the correct speed v,; at 
the front 7,(t) and followed by updating U(t, x,y). The new position of the front is 
stored implicitly in ¢;. With this approach, we avoid the difficulties that arise from 
explicitly tracking the front and thus increase the efficiency to deal with complex 
interfacial geometries. 

Given the normal speed, v;, at which the front 7,(t) moves, we would construct 


a speed function, Vi(t,2,y), which is a continuous extension of |vi(t, x, y)| from the 
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front 7,(t) over the whole computational domain. The governing equation of ¢, is 
then given by 
Ooi 


OE + Vi|Voi| = 0. (3.30) 


This equation will move ¢; with the correct speed at the front by assuring that 7,(t) 
will always coincide with the zero level set of ¢; at time t. 


We use ¢, to define the outward normal vector n, corresponding to 7, by 


From Equation (3.26) and (3.31), we can rewrite the expression for the speed of 7; (t) 


as 
V1 
|\Vé1| 


Since V, is equal to |v; (t)| along the interface, we can combine Equation (3.30) and 


vi(t) = —,|VU ny = —,|VU] (3.32) 


Equation (3.32) to get the following equation, which of course is only valid on the 
zero level set of @: 
Oo, 


Ge = MAIVUI|Voil, (wy) € 710) (3.33) 


Next, we need to extend the velocity function V; to a neighborhood of 7;(t). 
Step 2. Compute the velocity filed V(t, x, y) of U. 


In the algorithm, we compute approximations to Vi(t,z,y) at every grid point. 
One issue in computing VU arises from the fact that its approximation is usually in 
the order O(1) at points close to or on the front. By introducing V, defined as an 
extension of |vi(t,z, y)| = up| VU (Et, x, y)| from (x, y) € ™(t), we can avoid unnecessary 
numerical difficulties when we solve Equation (3.32) and (3.33). 

The approximation to VU at 7;(t) is based upon approximations to the derivatives 


of U in four coordinate directions to cut down on grid orientation effects. 
1 au\* (au? au\? (au? 
ronan $((B) SAB) 
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Each approximation to a derivative of U can be continuously extended away from 


the front by the advection equations 


Y 


Figure 3.5: Four coordinate directions used to compute the normal velocity. 


u; + S(t = 0, (3.34) 
ue+ S(615 2) =} (3.35) 
ue + (6158 = 0, (3.36) 
us + (ul = 0, (3.37) 


where u! = 0U/0z, u? = OU/Ody, u? = OU/On and u* = OU/OC¢ on 7(t). Here S is 
equal to the sign function. Equation (3.34) through Equation (3.37) have the effect of 
continuously extending u!, u?, u, u* away from the front by advecting these fields in 
the proper upwind direction. These equations will not degrade the value of V; on the 


front because ¢) is zero on 7;(t), so are S(¢ Gor) S(d1 ea S(d1 or) and S(¢1 sae 
Step 3. Update ¢, to be a signed distance function for one time step. 


From Equation (3.30) and (3.31), it is clear that the computation of the normal 
velocity, and normal vector at the front are all dependent upon the level set function 
@,. However, by Equation (3.30), level set function will cease to be an exact distance 


function even after one time step. In order to keep the accuracy of n,, and V;, we 
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need to avoid having steep or flat gradients developed in @,. One way to avoid these 
numerical difficulties is to reinitialize @,; to be an exact distance function from the 
evolving front 7,(t) at each time step. 

Firstly, we update ¢, by solving (3.30) with V, given in Step 2 and initial value 
¢1 obtained in Step 1. We then follow the same ideas as in [56] to reinitialize ¢y. 
In that paper, an algorithm was presented for reinitializing the level set function ¢ 
to be an exact signed distance function from the front. The basic idea behind this 
method is that given a function ¢? that is not a distance function, one can evolve it 
into a function ¢, that is an exact signed distance function from the zero level set of 


¢%. This can be accomplished by iterating the following equation to a steady state 


do 

Sr = S(¢1)G - Veal), (3.38) 
where ¢1(0,2,y) = ¢? and S again denotes the sign function. As in [56], the sign 
function S is smoothed by the equation 

oo 


Sg) a 
(1) iy pe 


(3.39) 


to avoid numerical difficulties while implemented. 

By using this approach, we avoid having to explicitly find the contour ¢? = 0 and 
then resetting values of the front ¢9 at grid points. From Equation (3.38), it is clear 
that the original position of the front will not change, but at points away from 7;(t), 


¢1 will be evolved into a distance function. 
Step 4. Update U(t, xz, y). 


After moving ¢, by the correct velocity at the front and then reinitializing ¢, to 
be an exact signed distance function from 7,(t) in Step 3, next we update U(t, x, y). 
Updating U(t, x,y) essentially boils down to solving the nonlinear parabolic partial 
difference equation (3.23) over the whole computational domain in the following four 


cases: 
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iil i+fj-1 


Figure 3.6: Illustration of updating U and V. 


e At points away from the front, which means the nearby four grid points are 
all inside the domain {2\(t), we solve the nonlinear parabolic partial difference 
equation (3.28) using a standard five-point stencil. For example, the grid point 


(i+ 1,7) in Figure 3.6. 


e For points near the front 7,(t), some special care should be taken. We effectively 
capture the front using the level set function @,. We can use one-sided different 
sign of ¢; to incorporate the distances between a point on the front and grid 
points neighboring it in either the vertical or horizontal direction. For example, 
yr © T(t), yp = (@— Dh, yr) € T1(t), we consider two grid points (i, 7 + 1) and 
(t,7) which border y,;. In y-direction, we have y; < yr < yj+i. We introduce 


or(i, J) 
bi(i, J “ir 1) -- br(i, J) 


jp=y STS 


and use U;,;, Uj;j;~-1, Uij—2, r and U(yr) = 0 to construct interpolating polyno- 


mial P. When we update U; ;, we use a standard five-point stencil by employing 


J? 
P(jh) instead of Ujj;41. For the case when front interacts with x-axis, we use 


the same process in x-direction. 


e For the extreme configuration, where we cannot find enough grid points inside 
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the domain to construct interpolating polynomial P, we employ the nearby grid 
points in the domain, and also the intersect points of the front and x and y-axis 


instead, to do Taylor expansion at that grid point to update U. 


e If a grid point lies on the front, we set the value U = 0 at that point (in view 


of (3.28)). For example, we set U;_1,;=0 for the grid point (i — 1,7). 


Step 5. Repeat Step 2 through Step 4 to @2 and V. 


Step 6. Repeat Step 2 through Step 6 to update ¢), ¢2, U and V for next time 


step. 


3.5 NUMERICAL EXPERIMENTS 


3.5.1 NUMERICAL TESTS FOR FRONT-TRACKING METHOD AND FRONT-FIXING METHOD 


OF 1D TWO-SPECIES COMPETITION-DIFFUSION MODEL 


Convergence test for front-tracking method of 1D competition-diffusion model 

In the 1D two-species competition-diffusion model (3.1)-(3.7) with parameters 
values (Dy, 1,71, K1, S°) =(0.4, 5, 2,1,0.4), (Do, ue, Yo, Ke, S9)=(0.4, 10, 1,2,1), L= 
1.2, Up = 2cos(*) and Vo = 4cos(%). Here we test the order of convergence in space 
with very refined temporal step size. 

In Tables 3.1-3.2, the error (both Lz and L,,) and the convergence to the solution 
of front-tracking method is examined, with final time t.,g = 0.01. The error is 
computed by the difference of the numerical solution with the exact solution. For 
all the examples below when the exact solution is not given, the solution with a 
fine resolution will be considered as reference or "exact" solution. As expected, a 
second-order convergence in space for both u and v can be observed. 


Convergence test for front-fixing method of 1D competition-diffusion model 
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Table 3.1: Accuracy results of species U for front-tracking method of 1D two-species 


model 


Table 3.2: Accuracy results of species V for front-tracking method of 1D two-species 


model 


N, xX N; || LoFrror | Order | L,.Error Order 
Accuracy test of U for front-tracking method 

61xle5 || 5.637e-04 1.699e-03 

121x1e5 || 1.035e-04 | 2.45 | 3.260e-04 2:88 
241x1e5 |) 1.850e-05 | 2.48 | 6.019e-05 2.44 

481x1e5 || 2.987e-06 | 2.63 | 9.833e-06 2.61 
961x1e5 || Reference 

Accuracy test of the front of U for front-tracking method 
61xled5 || 1.222¢e-04 2.233e-03 

121x1e5 |} 2.280e-05 | 2.42 | 5.672¢e-04 1.98 
241x1e5 |) 4.300e-06 | 2.39 | 1.296e-04 2.13 
481x1e5 || 8.000e-07 | 2.50 | 2.494e-05 2.38 
961x1e5 || Reference 


N, x Nz || LoError | Order | L,.Error Order 
Accuracy test of V for front-tracking method 

61x1le5 || 4.443e-04 1.373e¢-03 

121x1le5 || 7.882e-05 | 2.49 | 2.493e-04 2.46 

241x1e5 |) 1.396e-05 | 2.50 | 4.254e-05 2.55 

481x1e5 || 2.378e-06 | 2.55 | 9.871e-06 a 

961x1e5 | Reference 

Accuracy test of the front of V for front-tracking method 

61xle5 || 1.385e-04 3.268e-03 

121x1e5 || 2.721e-05 | 2.35 | 8.504e-04 1.94 

241x1e5 || 6.000e-06 | 2.19 | 1.922e-04 2.15 

481x1e5 || 1.200e-06 | 2.30 | 3.788e-05 2.34 

961x1e5 || Reference 


In the 1D two-species competition-diffusion model (3.1)-(3.7) with parameters val- 
ues (Di, 11,71; Ky, or) =(0.4, 5, 2. 1, 0.4), (Do, [2, 72; Ko, SSCA: 10, 1, 2, 1), Up = 
2cos(%) and Vo = 4cos(5). Here we test the order of convergence in space with very 


refined temporal step size. 
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In Tables 3.3-3.4, the error (both Zz and L,,) and the convergence to the solution 


of front-fixing method is examined, with final time tenqg = 0.01. 


As expected, a 


second-order convergence in space for both u and v can be observed. 


Table 3.3: Accuracy results of species U for front-fixing method of 1D two-species 


model 


Nz x N: L2Error | Order | L.Error Order 
Accuracy test of U for front-fixing method 

101x1e6 |) 1.195e-04 2.119e-04 

201x1e6 || 3.142e-05 | 1.93 | 5.424e-05 1.97 
401x1e6 |) 8.233e-06 | 1.93 | 1.3814e-05 2.05 
801x1e6 || 1.956e-06 | 2.07 | 3.983e-06 We 
1601 x1e6 || Reference 
Accuracy test of the front of U for front-fixing method 
101x1e6 || 3.178e-06 9.366e-05 

201x1e6 || 7.880e-07 | 2.01 | 2.424e-05 1.95 
401x1e6 |) 1.880e-07 | 2.07 | 5.927e-06 2.03 
801x1e6 || 3.800e-08 | 2.32 | 1.202c-06 2.30 
1601 x1e6 || Reference 


Table 3.4: Accuracy results of species V for front-fixing method of 1D two-species 


model 


Nz xX Nz L,Error | Order | L.Error Order 
Accuracy test of V for front-fixing method 

101xle6 | 1.038e-03 2.115e-03 

201x1e6 |) 2.776e-04 | 1.90 | 5.609e-04 1.91 
401x1le6 || 6.861e-05 | 2.02 | 1.381e-04 ZZ 
801x1le6 || 1.398e-05 | 2.30 | 2.807e-05 2.30 
1601x1e6 |) Reference 
Accuracy test of the front of V for front-fixing method 
101x1e6 |) 2.890e-05 7.595e-04 

201x1e6 |) 7.700e-06 | 1.91 | 2.123e-04 1.84 
401x1le6 || 1.910e-06 | 2.01 | 5.454e-05 1.96 
801x1e6 || 3.900e-07 | 2.29 | 1.141e-05 2.26 
1601x1e6 || Reference 


Comparison between front-fixing and front-tracking for 1D competition-diffusion model 
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In Figures 3.7-3.8, we use front-fixing method and front-tracking method to simu- 
late the 1D two-species competition-diffusion model (3.1)-(3.7) with parameters val- 
ues (D4, 41,71, K1, $9) =(0.4,5,2,1,0.4) , (Do, fe, ye, Ko, S2)= (0.4, 10,1,2,1), Uo = 
2cos(4"), and Vo = 4cos(%*) and spatial size h = 0.00125. It shows that the results 
of front-tracking method and the results of front-fixing method are consistent with 


each other. 


u(x,t) of front-tracking method V.S. u(x,t) with front-fixing method uH(t) of front-tracking method V.S. uH(t) of front-fixing method 
1.8 T T T T 0.5 T T 


- - — front-tracking - - — front-tracking 
front-fixing 0.49 front-fixing 


0.01 


Figure 3.7: Comparison numerical results of density of population and moving bound- 
ary of species U between front-tracking method and front-fixing method for 1D model 


v(x,t) with front-tracking method V.S. v(x,t) with front-fixing method vH(t) with front-tracking method V.S. vH(t) with front-fixing method 
4 T T T 7 1.14 : T T 
= = = front-tracking = = = front-tracking 


front-fixing 


front-fixing 


2.5/7 


=0.01) 
» 


Vv(x,t: 


Figure 3.8: Comparison numerical results of density of population and moving bound- 
ary of species V between front-tracking method and front-fixing method for 1D model 
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3.5.2. NUMERICAL DICHOTOMY: SPREADING VERSUS VANISHING 


Example 1. In 1D two-species competition-diffusion model (3.1)-(3.7) with pa- 
rameters values (D1, 41,71, K1, S9)=(2, 10, 2, 0.5, 4), (Da, fe, Ye, Ke, $9)=(1, 1, 1, 2, 1), 


Uo = cos(xe0) and Vo = 2cos(ze0): 


Figure 3.9 shows the case where the species U spreads successfully and V vanishes 


eventually. 
; The density of u at t=12 Z The free boundary of u 
09 22 
08+ 20 
O7- 18 
os, 16 
Mos eu 
ey << 
= 0.4, 12 
03, 10 
02, 8 
O1- 6 
0 L n 4 L L L L L 
0 5 10 15 20 25 0 2 4 6 8 10 12 
x t 
9x0" The density of v at t=12 i The free boundary of v 


Figure 3.9: Numerical simulation: species U spreads successfully and V vanishes 
eventually. 


Example 2. In this example, we take parameters values (Dj, (41,71, Ki, $?)= 
(4,2, 2,0.5,1), (Da, fe, ye, Ko, $2)= (2,1, 1,2, 0.6), Uo = cos( x0) and Vo = 3008( 300). 
Figure 3.10 shows the case where both species U and V vanish eventually. 

Example 3. In 1D two-species competition-diffusion model (3.1)-(3.7), parame- 


ters values are set as (Dj, 1,71, Ki, S?)=(0.5, 0.5, 1, 0.25, 2), (Do, fe, Yo, Ko, S3)= 


(1,10, 1,4, 2), Uo = 3008( 5am) and Vo = 2cos(xe0): Figure 3.11 shows the case where 
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The free boundary of u 


The density of u at t=10 


The free boundary of v 


The density of v at t=10 


Figure 3.10: Numerical simulation: both species U and V vanish eventually. 


both species U and V spread eventually. 


3.5.3 NUMERICAL TESTS FOR LEVEL SET METHOD FOR 2D MODEL WITH DIFFERENT 


INITIAL CONFIGURATION 


Example 4. In the 2D two species competition-diffusion model (3.23)-(3.29) with 
parameters values(D,, t1, 71, K1)=(4, 10, 1, 0.6), (Do, He, ¥2, K2)= (0.4, 5,3,0.5), the 
initial boundary of species U is set to be an equilateral triangle which centers at the 
origin point (0,0) with side length 1, while the initial boundary of species V is set 
to be a circle which centers at the origin point (0,0) with radius=1.5. The initial 
values Up(x,y), Vo(x, y) and the initial level set functions ¢9(x, y), d3(x, y) are set as 


following 
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The density of u at t=30 The free boundary of u 


0 5 10 15 20 25 30 36 40 0 5 10 15 20 25 30 


x t 


The density of v at t=30 The free boundary of v 


40(48 — 5 +) (v3 — y + Jg)(—V3x —y + 3g), (wy) € %(0), 


Uo(x, y) = 
0 (x,y) € Q5(0). 
A5cos(Y=F""), (a,y) € 22(0), 
Vo(z,y) = 
0 (x,y) € 25(0) 
V3a—-y+ V3ae—-y+4 
—min(B — 4, 4 y, Ae Ce (a,y) € (0), 
d(a,y)=4 0 (x,y) € O22,(0), 
. |\V3a—yt+--| |-v8a—-y+-+| sc 
min(|¥8 — 4, 4 yl, at EY (2,y) € M0) 


b5(a, y) = —(1.5 — y/2? + y?). 


Figure 3.12 shows the simulation of the evolvement of two species and their moving 


boundaries along time with an equilateral triangle as the initial boundary of U and 
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a circle as the initial boundary of V. In the figure of boundary line, the red curves 
represent the initial boundaries, and the blue curves simulate the evolvement of free 
boundaries. From Figure 3.12, we can see that the triangle evolves into a circle during 
the simulation. 

Example 5. In the 2D two species competition-diffusion model (3.23)-(3.29) with 
parameters values(Dj, t1, 71, 1)=(4, 20, 1, 0.6), (Da, fe, 2, K2)=(1, 5, 2, 0.5), the ini- 
tial boundary of species U is set to be is a square with side length=1, centered at 
(0,0), while the initial boundary of species V is set to be a circle which centers at the 
origin point (0,0) with radius=2. The initial values U(x, y), Vo(x, y) and the initial 


level set functions 9(x, y), 6§(z, y) are set as following 


101—a)1+a)1—y)+y), (ay) € 21(0), 


U(x, y) _ 
0 (x,y) € M§(0). 
eae 50cos(~——-—),_ (x, y) € 22(0), 
0 (x,y) € 23(0). 
—min(1l —|z|,1—lyl), (x,y) € 21(0), 
oi(x,y) = 4 0 (x,y) € 82,(0), 


min(|1— al], |1—Iyll) (ey) € 200), 


p5(x, y) = —(2 — x? + y?). 


Figure 3.13 shows the simulation of the evolvement of two species and their mov- 
ing boundaries along time with a square as the initial boundary of U and a circle as 
the initial boundary of V. In the figure of boundary line, the red curves represent 
the initial boundaries, and the blue curves simulate the evolvement of free bound- 
aries. From Figure 3.13, we can see that the square evolves into a circle during the 
simulation. 

Example 6. In the 2D two species competition-diffusion model (3.23)-(3.29) 
with parameters values (Dj, 1,71, K1)=(1, 15,1,0.4), (Do, Me, 2, K2)=(2, 5, 1, 0.5), 
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the initial boundary of species U is set to be is a circle which centers at (0,0) with 
radius=2.5, while the initial boundary of species V is set to be a circle which centers 
at the origin point (0,0) with radius=3.2. The initial values Uo(x,y), Vo(a,y) and 


the initial level set functions 6?(z, y), ¢8(x, y) are set as following 


40cos( ze (x,y) € 2,(0), 


Uo(z, y) -_ : 
0 (x,y) € £25(0). 
20cos(*¥ cay (x,y) € 22(0), 
Vo(x,y _ 
0 (x,y) € §25(0). 


(a, y) = —(2.5 — fa? + y?). 
08,9) ==B2= 2? +42). 


Figure 3.14 shows the simulation of the evolvement of two species and their moving 
boundaries along time with circles as the initial boundary of U and V. In the figure 
of boundary line, the red curves represent the initial boundaries, and the blue curves 
simulate the evolvement of free boundaries. From Figure 3.14, we can see that the 


circles propagate as circles during the simulation. 
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the contour line of u at =0.01 the contour line of v at =0.01 
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the contour line of v at t=0.04 
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Figure 3.12: The simulated dynamics where initial boundary of U is an equilateral 
triangle and initial boundary of V is a circle. The snapshots are taken at the times 
t = 0, 0.01, 0.04, respectively. 
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Figure 3.13: The simulated dynamics where initial boundary of U is a square and 
initial boundary of V is a circle. The snapshots are taken at the times t = 0, 0.05, 0.1, 


respectively. 
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Figure 3.14: The simulated dynamics where initial boundaries of U and V are circles. 


The snapshots are taken at the times t = 0,0.01, 0.05, respectively. 
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4.1 ABSTRACT 


In this chapter, we introduce Krylov implicit integration method to deal with the 
stiffness of the diffusive logistic model with a free boundary. We first transform 
the stiff reaction-diffusion system into a system with a fixed computational domain, 
and then introduce four different temporal schemes: Runge-Kutta, Crank-Nicolson, 
implicit integration factor (IIF) and Krylov IIF for handling such stiff systems. Nu- 
merical examples are examined to illustrate the efficiency, accuracy and consistency 
for different approaches, and it can be shown that Krylov IIF is superior to other 


three approaches in terms of stability and efficiency by direct comparison. 


4.2. INTRODUCTION 


According to the seminal paper by Du and Lin [14], the diffusive logistic model with a 
free boundary for the density of population of the invasive species U(t, x) depending 


on time ¢ and spatial variable x in one dimension can be described as follows: 


2 
= a = FW), t>0, 0<2< H(t), (4.1) 


together with the boundary conditions 


a (.0) =0, Ub) =O teal, (4.2) 


where the Stefan condition 


OU 
H'() = w(t. HD), > 0, (4.3) 
and the initial conditions 


The initial function U(x) satisfies the following properties: 
Uo(x) E C?((0, Hol), U;(0) = Uo( Ao) = 0, Uo(x) > 0, 0<a< Ho. (4.5) 
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Here H(t) is the unknown moving boundary such that the population is distributed in 
the interval |0, H(t)], and D > 0 is the dispersal rate. The parameter ps > 0 involved 
in the Stefan condition (4.3) is the proportionally constant between the population 
gradient at the front and the speed of the moving boundary. 

Few efficient numerical methods are proposed to solve the diffusive logistic model 
with a free boundary when the system is severely stiff. The difficulty is that extremely 
small time steps are required due to the stiffness of the system. When the explicit 
schemes are applied to solve such a system, due to stability constraints, an extremely 
small time step should be used and it might take a long time to finish one single 
simulation. However, while applying an implicit scheme [5, 25, 40] may be able to 
remove the stability constraints on the time step At, it usually requires solving a large 
global system of nonlinear equations for each time step, and the computational cost 
could be significant. To efficiently solve stiff reaction-diffusion equations, the authors 
in [43] have developed a class of implicit integration factor (IF) methods that are 
computationally much cheaper than fully implicit schemes and can be unconditionally 
stable or have generous stability conditions. The flexibility of representation of IIF 
method allows either direct calculation of the exponentials of matrices, or the use of 
Krylov subspace [10, 29, 31, 32, 39, 50] for non-constant diffusion coefficients or/with 
moving boundaries to compute their exponential matrix-vector multiplications for 
further saving in storage and computational cost. 

While very little existing work accounts for solving a stiff system with a free 
boundary, this chapter aims to develop and compare different numerical schemes to 
solve the system (4.1)-(4.5) in one dimension. We first use the Landau transforma- 
tion {12, 34] to convert the problem (4.1)-(4.5) into a system with fixed computational 
domain, in which the moving domain is included as a new separate variable coupled 
in the system. The Crank-Nicolson scheme and the Runge-Kutta method are applied 


to preserve second order accuracy both in time and space for better stability. We 
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extend the implicit integration factor (IIF) method to a system with a free boundary, 
which successfully removes the stability constraint associated with the diffusion and 
the stiff reaction terms thus allowing for large time step size. In order to calculate 
the matrix exponentials more efficiently, Krylov subspace incorporated with [IF is 
employed to speed up the simulation while maintaining the similar accuracy and sta- 
bility conditions. Numerical examples are performed to show the accuracy, stability 


and efficiency of four different proposed approaches. 


4.3 NUMERICAL SCHEMES FOR STIFF SYSTEMS WITH A FREE BOUNDARY 


As discussed in [48], the system (4.1)-(4.5) in one dimension can be transforming 
into a problem with a fixed computational domain [0,1]. Let us consider using the 


Landau transformation [12, 34], 


Zo) = WE, 2) = Ub, @). (4.6) 


H(t)’ 


Then the 1D diffusive logistic model (4.1) turns into the form, 


2 
OW ong ZW _ yp PW _ 


ey) Ot ( )5 Oz Oz? 


G(t)W(a—bW), t>0, O<z<1, (47) 


G(t) = H?(t), #20. (4.8) 


The boundary conditions (4.2) together with Stefan condition (4.3) take the following 


forms 
ar (60) =0, W(t,1)=0, t>0, (4.9) 
and 
Gi) = et Le. C0; (4.10) 


respectively, while the initial conditions of (4.4) become 
G(0) = Hj, W(0,z) =Wo(z) =Uo(zHo), O<2z<1. (4.11) 


76 


The smoothness of conditions of (4.5) for the initial function Up(x) are translated to 


Wo(z) as follows, 
Wo(z) € C7([0,1]), W (0) =Wo(1) =0, Wo(z)>0, O<z<1. (4.12) 


After the transformation, the new problem lies in solving the nonlinear system 
(4.7) in the fixed domain [0, 1] for the new variable z and unknown functions W,G. 
To numerically solve the system, let us consider the time step size k = At, and spatial 
mesh size h = Az = 1/M, where M is a positive integer denoting the total number 
of intervals in [0,1]. The mesh points (t”, z;) are denoted with t” = kn, n > 0, z; = 
jh, 0<j <M. Let us denote w? as the numerical approximation of W(t”, z;) at the 


mesh point (¢”, z;), ie., w? & W(t", z;), and let g” be the numerical value of G(t"). 


4.3.1 CRANK-NICOLSON SCHEME FOR IMPLICIT TEMPORAL DISCRETIZATION 


In this section, Crank-Nicolson scheme is applied to update G and W for each time 
step with the transformed system (4.7)-(4.10), which includes the following three 


steps: 
Step 1. Evaluate G(t) at t”*1, 


The new transformed Stefan condition (4.10) is discretized using first order for- 
ward approximation for G’(t), and second order approximation of ow (t, 1) with three 


points backward, 


oe = —* (3why —4wi, ,+why2), n>. (4.13) 
So G(t"*") can be evaluated by the following, 
GS GP he. 0. (4.14) 


Specially, we use g' = g° + kg'(0) to evaluate the starting value g!. 


Step 2. Update W(t, z) at t"*! by using Crank-Nicolson scheme. 


ws 


Let us consider the forward approximation of the time derivative of W, 


wrtt_ w= aw 
i eae (Eee: (4.15) 


and the second-order central approximation of each spatial derivative of W, 


WW AR no n n 2 
Wry — Wy_1 re OW Wry 2w; + Wy44 O-W 


eee d ~~ Pe: ; 
oh ae Oz (t 25); h2 O22 (t 25) (4 16) 


Based on the transformed Stefan condition (4.10) and from (4.13)-(4.14), G’(#) can 


be evaluated by the following, 
gt we —E (swt? — dwt, + wrt),, (4.17) 


Putting (4.15)-(4.17) all together into equation (4.7), we obtain 


k 252. Qh gn gnh2 
1 a wet — wet git ) puict — Qwett + wy ) 
99 oh g” + kg!” gt h2 
1 
+ 5 (wi (a DU; )ier wie? (a — bw?*")). (4.18) 


wheren>0, O0<j< M-—1. 
For j = 0, Equation (4.18) involves the fictitious value w”, at the ghost point 
(t”,—h). The value w”, can be evaluated from the discretization of the boundary 


condition (4.9), 
wi —w, 


7 


Such a nonlinear system can be solved by Newton iteration or other solvers. 


Step 3. Reevaluate G(t) at "7! using BDF2. 


A ee Qk 
a i ie yy, Owii —4uit + wht), n>. (4.19) 


We use BDF2 scheme to reevaluate G(t) at t”t! for n = 1, 2,3... in order to preserve 


the second order accuracy with O(k?) + O(h?). When n = 0, we use g' as a starting 
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value. The same process follows to select starting values. Remark Consistency and 


error analysis Let us consider the problems (4.6)-(4.12), and denote the Equations 
(4.7), (4.9) and (4.10) in the following form, 


aw @)zaW Dew : 
ENN a5 be Ug 


t>0, 0<2z<l, 


cw,G)= (1,0) =0, #>0, 


£3(W, G) = G'(t) + mr (t 1)=0, t>0, 


and after discretization, the equations are approximated by 


n+1 _ 


w Ue Deeg wag ND 4 = Dt oe as 
L id gh J I+1 j-1 1 D4 j It 
“Laut suet got auf +p 
99 2h grt ght th? 
= (wa bw"). (a = bw) S00, nS 0. 0S 7 = = 1, 
gs j j j 
La(w,g) =! =0, n>0, 
lee: eg 
L3(w, g) = i 7 (Sw —4wiy_, t whys) =9, n>0. 


Denoting the local truncation error T7(W, G) as 


T(1);(W, G) = L1(W7, G") — Li(W7, G"), 

T(2);(W, G) = L2(W;, G") — Lo(W7’, G"), 

T(3)j(W, G) = L3(W;', G") — £L3(W;, G"), 
Suppose we are given exact solutions W and G of problems (4.6)-(4.12) at t”, here 
Wr = W(t", z;) and G™ = Git”) 
Tr(W,G) = (T(1)?, T(2) 


gj? 


. According to [55], if the local truncation error 
”,7(3)) tend to zero as k + 0, h — 0, we say that the 


scheme L(w,g) = (Li(w,g), Lo(w, g), L3(w, g)) is consistent with problem L(W,G) = 
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(L1(W, G), Lo(W, G), £3(W, G)). With some calculation, we can get 


ko°w k? BW 
PO) (W,G) = > OP ere Ue 25) + ag (7 25) 

law 1oW 
sort 2) Ji = 3 at —(t es emer 

kew,., . Paw Low 

= Fae Ot Ege Ca) + a ape tH 

Law. ew. | Kaw 

rear a) tka et By + Hye (72 29) + Sa, 


J, and J» are defined as 


bn HE OW ne OU) STP OW og 
=~ 4G Ges OU Ga — Bam TD ot 11%) 
2; OW ph? ow 


ap ge on, n 4 
2 (t", 24) soe ex ("51 + OL). 


— OW nr | GUE) FG" (73) 
Ja =— Fay (Os 2) Geemsty GC) 

_ 4 OW n+l 2h? PW n+1 
Lose 3G (e) 023 (er 5a) 
2; OW nt 2uh? BW in k?G" (r: 

Hag i) sey aoe bse Ue) 
4 Oz SGPT 02 2G(t"+1) 
2, PBW a, . G(r) 
16 on 793) G(tr+1 

2 PW ong 5 HEM) GMs) 
46 08°? Ge) 260) 
D&W RG"(73) D hk? Ww. 

2 (t ae a) —~s n+1 4 (t +t 5) 

2 az 2G(t"+1) —- 2: 12G(#"+1) Az 
Df OW a PG) r ‘ 

25 eG on (t 3) 5a Gat + O(h*) + O(k*). 

“ ew, 
T(2)i (W, G) =, Az3 (t , 06) 
re k? BG ph? ew 
T(3)i (W, G) = = cy dt3 (Ta, = 3 23 (t , 67). 
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where the introduced parameters are given by 


t” <1, 72,73 <0", 

tr-} < TA < (or 

21 ~ 01,,09)03, 05S Za, 
ZM—2 < 64,67 <1. 


By 06S Bix 


By assuming that W(t, z) is fourth-order partial differentiable with respect to z and 
third-order differentiable with respect to t, and that G(t) is third-order differentiable, 


we can find that the local truncation error satisfies 
T(i)?(W, G) = O(K7) + O(h?), t= 1,2,3. 


4.3.2. RUNGE-KUTTA METHOD FOR EXPLICIT TEMPORAL DISCRETIZATION 


In this section, we introduce the Runge-Kutta scheme to solve the transformed system 
(4.7)-(4.12). For this approach, the evaluation of W,G from t” to t”*' consists of the 


following five steps: 
Step 1. Evaluate G(t"*!/?). 


The transformed Stefan condition (4.10) is discretized by using forward approxi- 
mation in time for G’(t) and three points backward approximation in space of owt, 1), 


1,653 
es —* (3why —4wi, 1+ whys), n> 0. (4.20) 
Thus G(t"*!/2) can be evaluated by the following, 


Kom 
ght? — gm + 59 , n>0. (4.21) 


Step 2. Update W(t, z) at t"+1/?. 
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Let us consider the forward approximation of all time derivatives, 


nt+1/2 on W n+1/2 n 
x age ——___~_ Git 4.22 


and the second-order central approximation in space, 


wi,-we, OW, Ue G2 gp OPW « 
cas ah : od Oz (t ane ! 72 J x Oz2 (t 124) (4.23) 


Combining (4.20)-(4.23) together, the equation (4.7) can be discretized by 


n+1/2 n n n mn, n n n 
5 — 5 Fj Wit T Wj-1 g' = poi = 207 + Why w"(a — bw”) 
k/2 2 2h og” ie Ca a j 


wheren>0, O0<j<M-—1. 


Step 3. Repeat Step 1 to evaluate G(t) at t’*! using the updated W and G at 
t"+1/2. Step 4. Evaluate W(t, z) at "+! using the Midpoint Method. 


m+1/2 _ a 6 _ Awrtl/? ca ney. aS 
n+l n n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 
we — wi Ree ne SOR et a ge wi + Wis 
k 9 Qgrtt/2h g gntl/2p2 
awh Pa buy), n>0, O<f<SM-1. (4.24) 
16%; 
n+1/2 n+1/2 
wnt! aye pA (MT Winn tay 
j gO" gnti/2\9 oh 
n+1/2 n+1/2 n+1/2 
Waal 2s + W541 n+1/2 n+1/2 
+ kD el ape t kw; "(a — bw," "*). 


Step 5. Reevaluate G(t) at t”*! using backward differentiation formula (BDF 2). 


g’t) = <g" — <9" 1 — =~ (3wip! —4uit, + wtf), n>. (4.25) 


remark Similar to Crank-Nicolson scheme, we have also performed consistency and 
error analysis to the Runge-Kutta scheme, which exhibits the second order accuracy 
in both space and time O(k?) + O(h?). We omit the details of derivation here to keep 


the paper more focused and concise. remark 
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4.3.3 IMPLICIT INTEGRATION FACTOR (IIF) METHOD 
In this section, the implicit integration factor (IIF) scheme is extended to solve the 
transformed Stefan problem (4.7)-(4.10). According to the equation (4.7), 


! 2 
ON ire yt, aw (ha) + Wit 2)(a—DWt,2)), eae 


(4.26) 
and notice that W(t, 1) = 0 from the boundary condition. 


First let us consider the central approximation of the spatial derivatives, 


W(t, 241) == Wt, zj-1) OW . 
os 2s 0<j<M-I1 
2h Oz (25); aa 


Wt, 25-4) ae 2W (t, Zi) + W(t, ea) aa OW 
he ~ Az? 


(t,z;), O<j<M-1. 


Putting this into the right hand side (RHS) of (4.26), and we obtain a semi-discretized 
ODE system 

dw —> —> 
Here W = (W(t, 20); W(t, 21); W(t, 22); ..W(t, z_1)), and C(t) is the approximation 


—> 
matrix for the spatial derivatives, R(W) is the nonlinear reaction term. For instance, 


C(t) at t = t” denoted by C,, with size M x M is given by 


—2D 2D 
grh2 grh? 0 0 0 

D__ zg” —~2D D_ | ag'” 0 
g@h2 Ahg” g?h2 gh? i 4hg” 

D__ zag'™ —2D D zo” 
0 g@h? 4hg” grh? grh2 + 4hg” 0 
0 . 0 
0 0 
0 ae: 0 D__ zm-19'" =~2D 
r Cos grh? A4hg” g@h? 


Again based on the transformed Stefan condition (4.10), g’” in the matrix C,, can be 


evaluated by the second order approximation of ow ct”, 1) with three points, 
Y (3w", — dw, + wt») & g”™ (4.28) 
i why — 4why-1 + Wi-2) © J : 
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Since C(t) is nonlinear, we can rewrite C(t) as 
Ct) = Cle") + Ele); (4.29) 


where €(t) is the correction term [31]. Thus we obtain 


dw —> — — 
Multiplying by the integration factor e~* on both sides, we get 
dw —> —> —> 
ere =e ™(C,W + E()W + R(W)). (4.31) 


—> — 
Let W,,4; denote as the numerical approximation for W(t"*’). After integrating it 


over one time step from t” to t”*! = t” + At, we have 


At 
Weng = ert, + elndt | T(E (E+ r)W(t” +7) + R(W(t" +7)))dr. (4.32) 
0) 


Following the same ideas as proposed in [43], to evaluate the integral, we approx- 

= = 
imate the nonlinear term e~7(E(t" + 7T)W(t” +7) + R(W(t” +7))) by an (r — 1)th 
order Lagrange polynomial p(r) with interpolation points at ¢”*', t",... #°'?7" for 


r-th order scheme. If we denote Fw +7))=E(t™+ r)wen +7)+ RW (t" +7)), 


we obtain 
as = Toe pee GANG 
pr) = > &“FWra) J] GAR SpE (4.33) 
tet j=l, jx J" 
Then we can get 
At —> a ay At t22 pA Gt 
—Cnt n iCn At 
eT EW (t" 47))dr So elke FW.) | ap 
[ ge 2 0 eth ag (j —i)At 


Putting all these together, and finally the r-th order IIF scheme is given by 


— — — = — 
Wr = Wy, + At(Qngi F (Wrst) + ¥> On-iF(Wn-i)) (4.34) 
i=0 
with the coefficients 
(@+1)CnAt pAt, T=? i 
Hie —_— | ee es ee ee) (4.35) 
aS 0 gala, jeg (FG — OA 
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More specially, the second order scheme (IIF2) is of the following form 


= <> — — 
Wri = Wr + At(anyiF(Wn41) + OnF(Wa)), (4.36) 
where 
1 1 
An = Caake An+1 = 5" 


—> 
After evaluating W at t”*', we reevaluate G(t) at t’*! using BDF 2 to preserve the 


second order accuracy for the front G, 


4 1 n— nm nr n 
gt = 37° ~ 39 mes 3, uM — dw + wife), n=l (4.37) 


remark The major computational cost of IIF arises from the evaluation of exponential 
matrices and multiplication of matrix with vectors. Since C,, is not a constant matrix, 


CuAt needs to be computed at every time step that is extremely 


its exponential e 
expensive thus will significantly slow down the simulation. Hence computation of 
exponential matrices at each time step should be avoided. Otherwise, it will become 


a bottleneck of regular IIF method for solving such system. remark 


4.3.4 KRYLOV SUBSPACE WITH IIF 


As discussed before [10], it is extremely expensive to construct matrix C,, and evaluate 


—> 
erMtW,, for each time step. Although the matrix C,, is sparse for most cases, the 


CnAt 


exponential matrix e is dense in general. To overcome this difficulty, applying 


AAt 


Krylov subspace for evaluation of the type of e*°'v is an excellent option. For 


example, in (20, 50], the Krylov subspace methods were used for the approximation of 
eA“ty, where A is a large sparse matrix and v is a given vector. Furthermore, Krylov 
subspace has been successfully incorporated with IIF and sparse grid structure for 
solving systems of PDEs with non-constant coefficient and unstructured grid meshes 
among many other systems {10, 31, 32, 39]. 

Following the literature (e.g. [20, 41]), we next illustrate how to apply the Krylov 


A 


subspace to evaluate e4*‘v in general. First the large sparse matrix A can be pro- 
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jected to the Krylov subspace 
Kg = span{v, Av, A’v,..., AC tv}. (4.38) 


The dimension Q of the Krylov subspace is usually much smaller than the dimension 
of the large sparse matrix A. For instance, we take Q = 25 for most of our test 
simulations in the next section. An orthonormal basis Vg = |v, V2, v3, ..., va] of the 
Krylov subspace Kg is generated by the well-known Arnoldi algorithm [57] as the 


following, 
1. Compute the initial vector: vy = v/||v||o. 
2. Perform iterations: Do 7 = 1,2,...,Q: 


a) Compute the vector w = Av,. 
b) Doi=1,2,...,j: 
i. Compute the inner product h;,; = (w, vi). 
ii. Compute the vector w = w — hy;v;. 
c) Compute hj41,; = ||w]|e. 
d) If hj41; = 0, then 
stop the iteration; 
else 


compute the next basis vector vj41 = w/hj41,;- 


In the Arnoldi algorithm, if h;,1,; = 0 for some 7 < Q, it means that the conver- 
gence has already occurred and the Krylov subspace is equal to span{vy, v2, ..., v;}- 
Hence the iteration can be stopped at this step 7, and we assign the value of this 7 
to Q. This algorithm will produce an orthonormal basis Vg of the Krylov subspace 
Kg. Denote the Q x Q upper Hessenberg matrix consisting of the coefficients h;; by 


Hg. Since the columns of Vg are orthogonal, we have 


Ha =VoAVo (4.39) 
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This means that the very small Hessenberg matrix Hg represents the projection of 
the large sparse matrix A to the Krylov subspace Kg, with respect to the basis Vo. 


Also since Vg is orthonormal, the vector VaVG eA“ty is the orthogonal projection of 


eA“'y on the Krylov subspace Kg, namely, it is the closest approximation to e44!v 


from Kg. Therefore 
eAAy we VaVge*'u = BVaVg 81 = BVQVg eA 'Vo er (4.40) 


where 6 = ||v|/2, and e; denotes the first column of the Q x Q identity matrix Jg. 


Noting from (4.39), we have the following approximation 
eAAty = BVgeZee, (4.41) 


Thus e4“*y for the matrix exponential with multiplication of a vector problem can 
be solved by a problem with much smaller size. The small matrix exponential e424* 
will be computed using a scaling and squaring algorithm with a Padé approximation 
with computational cost in the order of O(Q?), see [20, 26, 41]. 


Such ideas of Krylov subspace approximation can be easily applied to IIF scheme. 


For instance, the second order Krylov IIF scheme gives 


1 => = 
Wasa SAtF(Wngr) + BnWanetonAtey (4.42) 


= = s 
where 6, = ||W, + sAtF(W,,)|l2, Wan and Han are orthonormal basis and upper 
—> 
Hessenberg matrix generated by the Arnoldi algorithm with the initial vector W,, + 
—> 
sAtF(W,,). 


4.4 NUMERICAL EXPERIMENTS 


In this section, we will investigate the accuracy, stability and efficiency of the proposed 


four approaches: Runge-Kutta, Crank-Nicolson, IIF2 and Krylov IIF2, through a 
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number of test examples with the following form, 


ey _ pFU — u(a— dv), t>0, O0<a< H(t), 

2 (t,0) =0, U(t,H(t))=0, +t>0, (4.43) 
H'(t) = pS (t, A(t), t > 0, 

A(O) =p, 208) = oe); OS Alp: 


First we test the convergence order for all four methods in both space and time, 
in which the second order is observed for all cases. Next we show that the three 
schemes: Crank-Nicolson, IIF2 and Krylov IIF2 are much more stable than Runge- 
Kutta method. Finally Krylov IIF2 is shown to be able to dramatically reduce com- 


putational cost with comparison to Crank-Nicolson and IIF2. 


4.4.1 ACCURACY TEST 


Here we test the accuracy of the proposed Runge-Kutta scheme, Crank-Nicolson 


scheme, IIF2 scheme and Krylov IIF2 scheme in space and time by a simple case 


with D=a=b=p=1, Hy = 0.2, U0(x) = cos(3-) in the system (4.43). 


Since the exact solution is not known for this problem (4.43), for all our simu- 
lations, the numerical solution by using a very fine resolution will be considered as 
reference or "exact" solution. The differences between the numerical solutions and 
the "exact" solutions at the final time 7’ are measured by both Ly and L. errors. 
Here we set the final time 7’ = 0.1. We also compare all the "exact solution" from 
four different schemes to make sure that they are consistent with each other. 

In Tables (4.1-4.4), we run the numerical experiments at six different spatial and 
temporal resolutions. As expected, we can clearly see the second-order convergence 
rate in both space and time for all four schemes: Runge-Kutta, Crank-Nicolson, IIF2 


and Krylov IIF2. 
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Table 4.1: Accuracy results for Runge-Kutta method 


Nx N, || LofFrror | Order | LyError | Order 


Accuracy test of W 

26 x 5e4 1.85e-04 1.32¢e-04 
51x led 4.62e-05 2.00 | 3.28e-05 | 2.01 
101x2e5 1.16e-05 1.99 | 8.22e-06 | 2.00 
201 x4e5 2.89e-06 2.01 | 2.04e-06 | 2.01 
401 x 8e5 6.39e-07 2.18 | 4.50e-07 | 2.18 


801 x 16e5 Reference 
Accuracy test of G 
26x 5e4 2.66¢-04 6.09e-06 


51x1led5 6.65e-05 2.00 | 1.52e-06 | 2.00 
101x2e5 1.68e-05 1.99 | 3.85e-07 | 1.98 
201 x 4e5 4.20e-06 2.00 | 9.65e-08 | 2.00 
401x8e5 9.38e-07 | 2.16 | 2.17e-08 | 2.15 
801 16e5 Reference 


Table 4.2: Accuracy results for Crank-Nicolson method 


N,x N, || LoFrror | Order | LyError | Order 


Accuracy test of W 

26x 5e4 1.85e-04 2.68e-04 
51x1led 4.65e-05 | 2.00 | 3.30e-05 | 1.99 
101 x2e5 1.18e¢-05 1.98 | 8.30e-06 | 1.99 
201 x4e5 2.92e-06 2.01 | 2.06e-06 | 2.01 
401 8e5 6.30¢-07 2.21 | 4.38e-07 | 2.23 


801 x 16e5 Reference 
Accuracy test of G 
26x 5e4 1.33e-04 6.13e-05 


51x1led5 6.72e-05 2.01 | 1.54e-05 | 1.99 
101x2e5 1.71e-05 1.98 | 3.92e-06 | 1.97 
201 x 4e5 4.30e-06 1.99 | 9.90e-07 | 1.99 
401x8e5 9.30e-07 | 2.21 | 2.15e-07 | 2.20 
801 x 16e5 Reference 


4.4.2 STABILITY TEST 


In this section, we test the stability of the proposed four methods: Runge-Kutta 


scheme, Crank-Nicolson scheme, IIF2 scheme and Krylov IIF2 scheme by the system 
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Table 4.3: Accuracy results for IIF2 method 


N,x N, || LoFrror | Order | LyError | Order 
Accuracy test of W 

26 x 5e4 1.820-04 1.31e-04 
51x1e5 4.5le-05 | 2.02 | 3.20e-05 | 2.03 
101 x 2e5 1.11e-05 | 2.02 | 7.84e-06 | 2.03 
201 4e5 2.65e-06 | 2.07 | 1.86e-06 | 2.07 
401 8e5 5.34e-07 | 2.31 | 3.73e-07 | 2.32 


801 x 16e5 Reference 
Accuracy test of G 
26 x 5e4 2.65e-04 6.07e-06 


51x1e5 6.58e-05 2.01 | 1.51e-06 | 2.01 
101x2e5 1.64e-05 2.00 | 3.78e-07 | 2.00 
201 x 4e5 4.00e-06 2.04 | 9.26e-08 | 2.03 
401x8e5 8.35e-07 | 2.26 | 1.94e-08 | 2.26 
801 x 16e5 Reference 


Table 4.4: Accuracy results for Krylov IIF2 method 


N,x N, || Lofrror | Order | LyError | Order 
Accuracy test of W 

26 x 5e4 1.82¢e-04 1.31e-04 
51x1ed 4.5le-05 | 2.02 | 3.20e-05 | 2.03 
101 x2e5 1.11e-05 | 2.02 | 7.84e-06 | 2.03 
201 4e5 2.65e-06 | 2.07 | 1.86¢e-06 | 2.07 
401 8e5 5.30e-07 | 2.32 | 3.72e-07 | 2.32 


801 16e5 Reference 
Accuracy test of G 
26 x 5e4 2.65e-04 6.07e-06 


51x1le5 6.58e-05 2.01 | 1.51e-06 | 2.01 
101x2e5 1.64e-05 2.00 | 3.78e-07 | 2.00 
201 x 4e5 4.00e-06 2.04 | 9.26e-08 | 2.03 
401x8e5 8.40e-07 | 2.25 | 1.94e-08 | 2.26 
801 16e5 Reference 


(4.43) with parameters D = 1, = 0.50 = 2,6 = 1, Ho = 1.2, Ula) = cos( 7). For 
all simulations here, we set the final simulation time to Teng = 0.2, and the spatial 


grid N, = 801 while varying time step size. The reference solution is calculated by 
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choosing N, = 801 and At = 1078, and the maximum error is measured between the 
numerical solution and reference solution. 

In Figures (4.1) and (4.2), we can see that Runge-Kutta only converges for very 
small time step size, it quickly blows up as At increases, while all other three schemes: 
Crank-Nicolson, HF and Krylov IF allow for moderate large time step size, which 
exhibit extremely nice stability conditions. 

In Table 4.5, we also further test the order of accuracy for Crank-Nicolson, IIF2 
and Krylov IIF2 with fixed spatial resolution N,=801 and six different time resolution. 
The spatial resolution N,=801 is chosen fine enough such that the error is dominated 
by the time step. We list the maximum error between the numerical solution and the 
reference solution, and the order of accuracy. Clearly, the orders of accuracy in time 


for Crank-Nicolson, HF2 and Krylov IIF2 are two. 


—=-Crank-Nicolson 


-* KrylovilF 
-> Runge-Kutta 
poo TE 


-o 


max error 


10° 


107 107 


10° 10° 10° 


Lppitl L a 
107 
Timestep 


Figure 4.1: Error of U as a function of time step size. 


4.4.3. EFFICIENCY TEST FOR STABLE SCHEMES 


Since Crank-Nicolson, HF2 and Krylov IIF2 exhibit similar stability conditions, we 


test the efficiency for these three schemes in this section. To do this, here we consider 
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107 10% 10° 107 10° 107 107 
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Figure 4.2: Error of H as a function of time step size. 


Table 4.5: Error and order of accuracy in time for three stable schemes: Crank- 
Nicolson, IIF2 and Krylov IIF2 


At Crank-Nicolson IIF2 Krylov IIF2 
Accuracy test of W 
Ly. error | Order | L, error | Order | L, error | Order 
8.0 x 107° | 1.54e-8 - 6.50e-8 - 6.50e-8 - 
4.0 x 107° | 4.09e-9 1.91 2.23e-8 1.a5 2.23e-8 1.55 
2.0 x 10-° | 1.05e-9 1.97 6.28e-9 1.83 6.28¢-9 1.83 
1.0 x 107° | 3.22e-10 | 1.70 1.30e-9 227 1.30e-9 2.27 
5.0 x 107 | 8.14e-11 | 1.98 | 2.85e-10 | 2.20 | 2.85e-10 | 2.20 
2.5 x 10-6 Reference Reference Reference 
Accuracy test of G 
Log error | Order | Loo error | Order | La error | Order 
8.0 x 107° | 4.16e-7 “ 8.49e-7 8.49e-7 Z 
4.0x107-° | 1.44e-7 1.53 2.81e-7 1.59 2.81e-7 1.59 
2.0x 107° | 4.86e-8 Lot 9.22e-8 1.61 9.22¢-8 1.61 
1.0x 107° | 1.54e-8 1.66 2.88e-8 1.68 2.88¢-8 1.68 
5.0 x 10-° | 3.92e-9 1.97 7.27e-9 1.99 7.27e-9 1.99 
25 10-° Reference Reference Reference 


two cases: case 1: large diffusion system with D = 100,a = b = 1; case 2: stiff system 
with D = 1,a = b = 100. For both cases, we set w = 1, Hp = 2, the final simulation 


time Teng = 0.1, and the initial condition Up = cos( 37). For all the simulations here, 
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we select three different spatial resolutions N, = 1001, N, = 2001 and N, = 4001 
with fixed time step size At = 107+. 


Table 4.6: CPU time for three stable schemes: Crank-Nicolson, HF2 and Krylov IIF2 
of large diffusion system. The unit of CPU time is second. 


At = 10-4 N, = 1001 N, = 2001 N, = 4001 

Crank-Nicolson 16.576 75.092 395.512 

Krylov IIF2 14.595 59.566 295.958 

IIF2 211.227 1099.695 9694.277 
25 x10 


-- Crank-Nicolson 
— KrylovilF 


0.5 | 


Figure 4.3: Numerical simulation: U and H for the large diffusion system. 


From Figures 4.3 and 4.4, first we show that the solutions agree with each other 
for three different schemes. We can also see that when compared to the stiff system, 
the front H for the large diffusion system moves much faster initially and quickly 
reaches the steady state afterwards. 

From Table (4.6), for the large diffusion system, it is clear that the regular IIF2 
is more than 20 times slower for the refined grid than the other two due to the very 
expensive evaluation of exponential matrices at every time step, while Krylov HIF 2 is 


slightly faster than Crank-Nicolson for all different grid resolution. 
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Table 4.7: CPU time for three stable schemes: Crank-Nicolson, HF2 and Krylov IIF2 
of stiff system. The unit of CPU time is second. 


fits i0-* N, = 1001 N, = 2001 N, = 4001 
Crank-Nicolson 30.149 141.601 760.832 
Krylov IIF2 16.706 71.501 346.278 
IIF2 389.514 1877.676 22626.109 


2.35 Crnak-Nicolson 
—KrylovilF 
23 IF 


H(t) 


Figure 4.4: Numerical simulation: U and H for the stiff system. 


From Table (4.7), we can show that all the simulations for the stiff system are 
much slower than those of large diffusion system with the same grid size. Once again, 
the regular IIF2 is extremely slower than the other two, while Krylov IIF2 is more 
than 2 times faster than Crank-Nicolson and such trend is more obvious when the 
grid becomes finer. To our understanding, for the Crank-Nicolson scheme, it requires 
solving nonlinear system for each time step, thus demanding more CPU time for the 
stiff systems. However, Krylov IIF2 avoids solving the nonlinear system, instead it 
only handles nonlinear equation for each grid point. Thus Krylov IIF can reduce 
computational cost for the system with very stiff reaction terms. In brief, applying 
Krylov subspace method incorporated with IIF is able to dramatically further improve 


the efficiency for solving very stiff systems with moving boundaries. 
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CHAPTER 5 


CONCLUSIONS 


In this dissertation, we are devoted to designing efficient numerical methods to in- 
vestigate the spreading-vanishing dichotomy of the diffusive logistic model with a 
free boundary and the spreading behavior of two invasive species modeled by a diffu- 
sive competition system with free boundaries. Combining the implicit solver, front- 
tracking method, which explicitly tracks the locations of interfaces, is adopted to 
develop numerical schemes for one dimension and higher dimensions with spherical 
symmetry. In particular, a front-fixing approach is also introduced to verify the accu- 
racy and stability of front-tracking method. In higher dimensional cases, we introduce 
level set method to handle topological bifurcations. It is demonstrated that level set 
method is able to robustly and efficiently capture different complicated geometries. 
Illustration with numerical examples of the spreading-vanishing dichotomic behavior 
and the spreading behavior of two invasive species are given. Numerical experiments 
in 1D and 2D spaces are also performed to demonstrate the accuracy, consistency 
and stability of the proposed numerical methods. 

We also derive four different numerical schemes: Runge-Kutta, Crank-Nicolson, 
IIF2 and Krylov IIF2, to systematically solve stiff reaction-diffusion system with 
a free boundary. Through numerical experiments, we show that the three schemes: 
Crank-Nicolson, IIF and Krylov IF exhibit very nice stability property, while Runge- 
Kutta only admits for very small time steps with the fine mesh size. Finally through 
efficiency test, it reveals that the regular IIF is too expensive for solving such systems 


due to the evaluation of exponential matrices at every time step, while Krylov IF is 


95 


the most efficient approach, especially for the system with very stiff reaction terms. 

Currently we are incorporating level set method with explicit exponential time 
differencing method for solving reaction-diffusion systems with moving boundaries in 
two dimensions, and some nice preliminary results have been obtained that further 
demonstrates the great promise of the proposed approach. Compared to 1D prob- 
lems, there are quite a few new numerical challenges when explicit exponential time 
differencing method is coupled with level set method for solving 2D problems. Firstly, 
due to the moving boundaries, the associated exponential matrices for discretized dif- 
fusion operator need to be assembled with caution by an ordered vector for each time 
step. Secondly, level set function should be carefully incorporated into explicit ex- 
ponential time differencing method solver for reaction-diffusion equations with some 
special numerical treatment, especially for the points near moving boundaries. Fi- 
nally, new ideas for extrapolation on the boundary points with different boundary 


conditions need to be explored to keep the higher order accuracy. 


96 


[10] 


BIBLIOGRAPHY 


D. G. Aronson and H. F. Weinberger, Nonlinear diffusion in population genetics, 
combustion, and nerve pulse propagation. In Partial differential equations and 
related topics, Springer, Berlin, Heidelberg, 1975, pp. 5-49. 


D. G. Aronson and H. F. Weinberger, Multidimensional nonlinear diffusion aris- 
ing in population genetics. Advances in Mathematics, 30 (1978), 33-76. 


W. Bao, Y. Du, Z. Lin and H. Zhu, Free boundary models for mosquito range 
movement driven by climate warming. Journal of Mathematical Biology, 76 
(2018), 841-875. 


G. Bunting, Y. Du and K. Krakowski, Spreading speed revisited: analysis of a 
free boundary model. Networks and Heterogeneous Media, 7 (2012), 583-603. 


K. Burrage and J. C. Butcher, Stability criteria for implicit Runge-Kutta meth- 
ods. SIAM Journal on Numerical Analysis, 16 (1979), 46-57. 


L. A. Caffarelli and S. Salsa, A Geometric Approach to Free Boundary Problems, 
American Mathematical Soc., 2005. 


Y. Cao, A. Faghri and W. S. Chang, A numerical analysis of Stefan problems 
for generalized multi-dimensional phase-change structures using the enthalpy 
transforming model, International Journal of Heat and Mass Transfer, 32 (1989), 
1289-1298. 


H. Chen, C. Min and F. Gibou, A numerical scheme for the Stefan problem on 
adaptive Cartesian grids with supralinear convergence rate, Journal of Compu- 
tational Physics,228 (2009), 5803-5818. 


S. Chen, B. Merriman, S. Osher and P. Smereka, A simple level set method for 
solving Stefan problems, Journal of Computational Physics, 135 (1997), 8-29. 


S. Chen and Y. Zhang, Krylov implicit integration factor methods for spatial 
discretization on high dimensional unstructured meshes: application to discon- 


97 


[15] 


[16] 


[17] 


[18] 


[19] 


[20] 


[21] 


tinuous Galerkin methods, Journal of Computational Physics, 230 (2011), 4336- 
4352. 


I. L. Chern, J. Glimm, O. McBryan, B. Plohr and S. Yaniv, Front tracking for 
gas dynamics, Journal of Computational Physics, 62 (1986), 83-110. 


J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984. 


Y. Du and Z. Guo, The Stefan problem for the Fisher-KPP equation, Journal 
of Differential Equations, 253 (2012), 996-1035. 


Y. Du and Z. Lin, Spreading-vanishing dichotomy in the diffusive logistic model 
with a free boundary, SIAM Journal on Mathematical Analysis, 42 (2010), 377- 
405. 


Y. Du and B. Lou, Spreading and vanishing in nonlinear diffusion problems 
with free boundaries, Journal of the European Mathematical Society, 17 (2015), 
2673-2724. 


Y. Du, H. Matano and K. Wang, Regularity and asymptotic behavior of nonlinear 
Stefan problems, Archive for Rational Mechanics and Analysis, 212 (2014), 957- 
1010. 


Y. Du and C.H. Wu, Spreading with two speeds and mass segregation in a 
diffusive competition system with free boundaries, Calculus of Variations and 
Partial Differential Equations, 57 (2018). https://doi.org/10.1007/s00526-018- 
1339-5. 


R. Fedkiw and S. Osher, Level Set Methods and Dynamic Implicit Surfaces, 
Applied Mathematical Sciences, 153. Springer-Verlag, New York, 2003. 


R. A. Fisher, The wave of advance of advantageous genes. Annals of eugenics, 7 
(1937), 355-369. 


E. Gallopoulos and Y. Saad, Efficient solution of parabolic equations by Krylov 
approximation methods, SIAM Journal on Scientific and Statistical Computing, 
13 (1992), 1236-1264. 


F. Gibou and R. Fedkiw, A fourth order accurate discretization for the Laplace 
and heat equations on arbitrary domains, with applications to the Stefan prob- 
lem, Journal of Computational Physics, 202 (2005), 577-601. 


98 


[22| 


[25] 


[26] 


(27 


uo 


[28] 


[29] 


[30] 


[31] 


[32] 


J. Glimm, X. L. Li, Y. Liu and N. Zhao, Conservative front tracking and level 
set algorithms, Proceedings of the National Academy of Sciences, 98 (2001), 
14198-14201. 


J.S. Guo and C.H. Wu, Dynamics for a two-species competition-diffusion model 
with two free boundaries, Nonlinearity 28 (2015), 1-27. 


J.S. Guo and C.H. Wu, On a free boundary problem for a two-species weak 
competition system, Journal of Dynamics and Differential Equations, 24 (2012), 
873-895. 


E. Hairer and G. Wanner, Stiff differential equations solved by Radau methods, 
Journal of Computational and Applied Mathematics, 111 (1999), 93-111. 


N. J. Higham, The scaling and squaring method for the matrix exponential 
revisited, SIAM Journal on Matrix Analysis and Applications, 26 (2005), 1179- 
1193. 


J. Hilditch and P. Colella, A front tracking method for compressible flames in 
one dimension, SIAM Journal on Scientific Computing, 16 (1995), 755-772. 


D. Hilhorst, M. lida, M. Mimura and H. Ninomiya, A competition-diffusion 
system approximation to the classical two-phase Stefan problem, Japan Journal 
of Industrial and Applied Mathematics, 18 (2001), 161-180. 


M. Hochbruck and C. Lubich, On Krylov subspace approximations to the matrix 
exponential operator, STAM Journal on Numerical Analysis, 34 (1997), 1911- 
1925. 


J. Hua, J. F. Stene and P. Lin, Numerical simulation of 3D bubbles rising in 
viscous liquids using a front tracking method, Journal of Computational Physics, 
227 (2008), 3358-3382. 


T. Jiang and Y. Zhang, Krylov implicit integration factor WENO methods for 
semi-linear and fully nonlinear advection-diffusion-reaction equations, Journal 
of Computational Physics, 253 (2013), 368-388. 


T. Jiang and Y. Zhang, Krylov single-step implicit integration factor WENO 
method for advection-diffusion-reaction equations, Journal of Computational 
Physics, 311 (2016), 22-44. 


99 


[33] 


[34] 


[35] 


[36] 


[37] 


[38] 


[39] 


[40] 


[41] 


[42] 


[43] 


[44] 


A. N. Kolmogorov, I. G. Petrovsky and N. S. Piskunov, Investigation of the equa- 
tion of diffusion combined with increasing of the substance and its application 
to a biology problem. Bull. Moscow State Univ. Ser. A: Math. Mech 1, 6 (1937), 
1-25. 


H. G. Landau, Heat conduction in a melting solid, Quarterly of Applied Mathe- 
matics, 8 (1950), 81-94. 


R. J. Leveque and Z. Li, The immersed interface method for elliptic equations 
with discontinuous coefficients and singular sources, SIAM Journal on Numerical 
Analysis, 31 (1994), 1019-1044. 


R. J. Leveque and K. M. Shyue, Two-dimensional front tracking based on high 
resolution wave propagation methods, Journal of Computational Physics, 123 
(1996), 354-368. 


X. Liang and X. Q. Zhao, Asymptotic speeds of spread and traveling waves for 
monotone semi-flows with applications, Communications on Pure and Applied 
Mathematics, 60 (2007), 1-40. 


S. Liu and X. Liu. Numerical methods for a two-species competition-diffusion 
model with free boundaries, Mathematics, 2018, 6, 72. 


D. Lu and Y. Zhang, Krylov integration factor method on sparse grids for high 
spatial dimension convection-diffusion equations, Journal of Scientific Comput- 
ing, 69 (2016), 736-763. 


M. M. Mac Low and R. S. Klessen, Control of star formation by supersonic 
turbulence, Reviews of Modern Physics, 2004, 76, 125. 


C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential 
of a matrix, twenty-five years later, SIAM Review, 45 (2003), 3-49. 


Q. Nie, F. Y. Wan, Y. Zhang and X. Liu, Compact integration factor methods 
in high spatial dimensions, Journal of Computational Physics, 277 (2008), 5238- 
5250. 


Q. Nie, Y. Zhang and R. Zhao, Efficient semi-implicit schemes for stiff systems, 
Journal of Computational Physics, 214 (2006), 521-537. 


S. Osher and R. P. Fedkiw, Level set methods: an overview and some recent 
results, Journal of Computational Physics, 169 (2001), 463-502. 


100 


[45] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: 
Algorithms based on Hamilton-Jacobi formulations, Journal of Computational 
Physics, 79 (1988), 12-49. 


[46] D. Peng, B. Merriman, S. Osher, H. Zhao and M. Kang, A PDE-based fast local 
level set method, Journal of Computational Physics, 155 (1999), 410-438. 


[47] C. S. Peskin, The immersed boundary method, Acta Numerica, 11 (2002), 479- 
O17. 


[48] M. A. Piqueras, R. Company and L. Jodar, A front-fixing numerical method 
for a free boundary nonlinear diffusion logistic population model, Journal of 
Computational and Applied Mathematics, 309 (2017), 473-481. 


[49] L. I. Rubinstein, The Stefan Problem, Providence, RI: American Mathematical 
Society, 1971. 


[50] Y. Saad, Analysis of some Krylov subspace approximations to the matrix expo- 
nential operator, SIAM Journal on Numerical Analysis, 29 (1992), 209-228. 


[51] J. A. Sethian, A fast marching level set method for monotonically advancing 
fronts, Proceedings of the National Academy of Sciences, 93 (1996), 1591-1595. 


[52] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving In- 
terfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and 
Materials Science, Cambridge University Press, 1999. 


[53] N. Shigesada and K. Kawasaki, Biologicallnvasions: Theory and Practice, in: 
Oxford Series in Ecology and Evolution, Oxford University Press, Oxford, 1997. 


[54] J. G. Skellam, Random dispersal in theoretical populations, Biometrika, 38 
(1951), 196-218. 


[55] G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Dif- 
ference Methods, Oxford University Press, 1985. 


[56] M. Sussman, P. Smereka and S. Osher, A level set approach for computing 
solutions to incompressible two-phase flow, Journal of Computational Physics, 
114 (1994), 146-159. 


[57| L. N. Trefethen and III D. Bau, Numerical Linear Algebra, SIAM, 1997. 


101 


[58] S. O. Unverdi and G. Tryggvason, A front-tracking method for viscous, incom- 
pressible, multi-fluid flows, Journal of Computational Physics, 100 (1992), 25-37. 


[59] M. Wang and J. Zhao, Free boundary problems for a Lotka-Volterra competition 
system, Journal of Dynamics and Differential Equations, 26 (2014), 655-672. 


[60] M. Wang and Y. Zhang, Note on a two-species competition-diffusion model with 
two free boundaries, Nonlinear Analysis, 159 (2017), 458-467. 


[61] H. F. Weinberger, On spreading speeds and traveling waves for growth and 
migration models in a periodic habitat, Journal of Mathematical Biology, 45 
(2002), 511-548. 


[62] H. F. Weinberger, M. A. Lewis and B. Li, Anomalous spreading speeds of coop- 
erative recursion systems, Journal of Mathematical Biology, 55 (2007), 207-222. 


[63] A. Wiegmann and K. P. Bube, The immersed interface method for nonlinear 
differential equations with discontinuous coefficients and singular sources, STAM 
Journal on Numerical Analysis, 35 (1998), 177-200. 


[64] C.H. Wu, The minimal habitat size for spreading in a weak competition system 
with two free boundaries, Journal of Differential Equations, 259 (2015), 873-897. 


[65] J. J. Xu, Z. Li, J. Lowengrub and H. Zhao, A level-set method for interfacial 
flows with surfactant, Journal of Computational Physics, 212 (2006), 590-616. 


[66] H. K. Zhao, T. Chan, B. Merriman and S. Osher, A variational level set approach 
to multiphase motion, Journal of Computational Physics, 127 (1996), 179-195. 


[67| L. Zhu and C. S. Peskin, Simulation of a flapping flexible filament in a flowing 
soap film by the immersed boundary method, Journal of Computational Physics, 
179 (2002), 452-468. 


102 


APPENDIX A 


PERMISSION TO REPRINT 


A.1 COPYRIGHT RELEASE BY INTERNATIONAL JOURNAL OF COMPUTER 


MATHEMATICS 
< E @ s100.copyright.com (6) (+) fo Oo 
Microsoft Office Home | Mail - LIU, SHUANG - Outlook Numerical studies of a class of reaction... Rightslink® by Copyright Clearance Cen... + 
Copyright . . ® wn 2 rl e e 
Clearance R | g hts Li al k A : e ea ~, 
< Center Home Help Live Chat Sign in Create Account 


Numerical studies of a class of reaction-diffusion equations with 
Stefan conditions 


Author: Shuang Liu, , Yihong Du, et al 
(w) Taylor & Francis publication: International Journal of Computer Mathematics 
ayer EeFranicis Grou Publisher: Taylor & Francis 
Date: Apr 4, 2019 


Rights managed by Taylor & Francis 


Thesis/Dissertation Reuse Request 


Taylor & Francis is pleased to offer reuses of its content for a thesis or dissertation free of charge contingent on 


resubmission of permission request if work is published. 
CLOSE 


| 


© 2019 Copyright - All Rights Reserved | Copyright Clearance Center, Inc. | Privacystatement | Terms and Conditions 
Comments? We would like to hear from you. E-mail us at customercare@copyright.com 


103 


A.2 COPYRIGHT RELEASE BY MATHEMATICS 


< f 


@ mdpi.com 


@ 


Thesis & Dissertation - Graduat.. 


https://www.sc.edu/study/colleg... 


Microsoft Office Home 


For Librarians 


For Publishers 


For Societies 


Article Processing Charges 


Open Access Policy 


Institutional O.A. Program 


Editorial Process 


MDP! | Information for Authors 


Favorites 


Manuscript Submission & Instructions for Authors 


Choose a journal... 


[ instructions for Authors ]| Aims & Scope I Submit Manuscript | 


Authors and Readers Benefit 
from MDPI’s Pledges to: 


publish thoroughly peer-reviewed journals of high 
scholarly impact 


maintain a quick publication procedure — manuscripts 
are published within 6-8 weeks of submission, provided 
that no major revisions are required 


publish full open access journals — readers can 
access all articles published on this platform for free, 
including state-of-the-art review articles 


publish citation-tracked journals — MDP! is working 
towards quick coverage and citation-tracking of all of its 
journals in the Science Citation Index Expanded (SCIE, 

see the list of covered journals) and Scopus (see the 

list of covered journals). In 2017, 79% of all published 
articles were covered by the SCIE, 


for biomed titles, MDP! is working towards quick 
coverage in PubMed Central and PubMed (see the list 
of covered journals). In 2017, 48% of all published 
articles were covered by PubMed. 


For Authors and Readers Open 
Access Means: 


© free availability of the literature on the Internet without 
any subscription or price barriers 


immediate open access once an article is released 
(no embargo period) 


« authors retain all copyrights - authors will not be 
forced to sign any copyright transfer agreements 


permission of re-use the published material if proper 
accreditation is given (Creative Commons Attribution 
License @) 


Read the full open access information here. 


Avoiding Surcharges from Extensive English Editing or Extensive 


104 


Back to Top 


A.3 


ie 


COPYRIGHT RELEASE BY DISCRETE AND CONTINUOUS DYNAMICAL 


SYSTEMS-SERIES B 


li @ outlook.office.com G ° fh ia) 


Thesis & Dissertation - Graduate School |... https://www.sc.edu/study/colleges_school... Mail - LIU, SHUANG - Outlook Format thesis workshop - Google 3¢#4 ar 


ch Print =X Cancel 


Ww’ be 

NG 
From: LIU, SHUANG [mailto:shuang@email.sc.edu] 

z. Sent: Wednesday, November 13, 2019 2:02 PM 
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" To: editorial@aimsciences.org 

3 Subject: Author's copyright release 

wo 
Dear Editorial Manager, 

ii] 
I'm the first author of Krylov implicit integration factor method for a class of stiff reaction-diffusion systems with free 
boundaries which is accepted in Discrete and Continuous Dynamical Systems-Series B. 

J https://www.aimsciences.org/article/doi/10.3934/dedsb.2019176 

Q I'm writing to get the permission to reuse part content of the article as part of my PhD dissertation. The work in the article is part of 
my research during my PhD program. 

2 
Please let me know if you need more information to process the permission. Thanks a lot for your help! 

> 
Best Regards, 

i Shuang 
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Cc: ‘editorial’ <editorial@aimsciences.org> 


Dear Shuang, 
You have the permission provided the proper acknowledgement is made of the original publication. Thanks. 


Best regards, 
Liwei Ning 


Liwei Ning 

Editorial Manager 

American Institute of Mathematical Sciences 
Tel: 417-351-3204 

Fax : 417-351-3204 

Email: editorial@ aimsciences.org 
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